An approximate empirical Bayesian method for large-scale linear-Gaussian inverse problems

TL;DR

This paper introduces an efficient approximate Bayesian approach for large-scale linear-Gaussian inverse problems, leveraging low-rank and randomized methods to evaluate marginal likelihoods, enabling practical hyperparameter estimation.

Contribution

It proposes a novel low-rank approximation technique for marginal likelihood evaluation, improving computational feasibility in large-scale Bayesian inverse problems.

Findings

Accurately approximates marginal likelihoods using low-rank methods

Demonstrates efficiency with randomized SVD and spectral approximations

Shows good performance in numerical experiments

Abstract

We study Bayesian inference methods for solving linear inverse problems, focusing on hierarchical formulations where the prior or the likelihood function depend on unspecified hyperparameters. In practice, these hyperparameters are often determined via an empirical Bayesian method that maximizes the marginal likelihood function, i.e., the probability density of the data conditional on the hyperparameters. Evaluating the marginal likelihood, however, is computationally challenging for large-scale problems. In this work, we present a method to approximately evaluate marginal likelihood functions, based on a low-rank approximation of the update from the prior covariance to the posterior covariance. We show that this approximation is optimal in a minimax sense. Moreover, we provide an efficient algorithm to implement the proposed method, based on a combination of the randomized SVD and a spectral approximation method to compute square roots of the prior covariance matrix. Several numerical examples demonstrate good performance of the proposed method.

Equations

$$\mathbf{y}=G\mathbf{x}+\mathbf{\eta }$$

$$\pi (\mathbf{x}|\mathbf{y})=\frac{\pi (\mathbf{y}|\mathbf{x})\pi (\mathbf{x})}{\pi (\mathbf{y})}.$$

$$\mathbf{y}\mid \mathbf{x}\sim \mathcal{N}(G\mathbf{x},{\Gamma }_{\mathrm{obs}}),\mathbf{x}\sim \mathcal{N}(0,{\Gamma }_{\mathrm{pr}}).$$

$$\mathbf{x}\mid \mathbf{y}\sim \mathcal{N}({\mathbf{\mu }}_{\mathrm{pos}},{\Gamma }_{\mathrm{pos}}),$$

$${\mathbf{\mu }}_{\mathrm{pos}}={\Gamma }_{\mathrm{pos}}{G}^{\top }{\Gamma }_{\mathrm{obs}}^{-1}\mathbf{y}\text{mbox}and{\Gamma }_{\mathrm{pos}}={\left(H+{\Gamma }_{\mathrm{pr}}^{-1}\right)}^{-1},$$

$$H=G^{\top }{\Gamma }_{\mathrm{obs}}^{-1}G$$

$$\underset{\mathbf{\theta }}{\max }\pi (\mathbf{y}|\mathbf{\theta })=\int \pi (\mathbf{y}|\mathbf{x},\mathbf{\theta })\pi (\mathbf{x}|\mathbf{\theta })d\mathbf{x}$$

$$\underset{\mathbf{\theta }}{\max }\pi (\mathbf{\theta }|\mathbf{y})=p(\mathbf{\theta })\int \pi (\mathbf{y}|\mathbf{x},\mathbf{\theta })\pi (\mathbf{x}|\mathbf{\theta })d\mathbf{x}$$

$$\underset{\mathbf{\theta }}{\min }-\log \pi (\mathbf{y}|\mathbf{\theta })$$

$$=\frac{1}{2}{\mathbf{y}}^T{\Gamma }_{\mathrm{obs}}^{-1}\mathbf{y}+\frac{1}{2}\log |{\Gamma }_{\mathrm{obs}}|-\frac{1}{2}{\mathbf{z}}^T{\Gamma }_{\mathrm{pos}}\mathbf{z}+\frac{1}{2}\log \frac{|{\Gamma }_{\mathrm{pr}}|}{|{\Gamma }_{\mathrm{pos}}|},$$

$$L(\mathbf{\theta },\mathbf{z}):=-\frac{1}{2}{\mathbf{z}}^T{\Gamma }_{\mathrm{pos}}\mathbf{z}+\frac{1}{2}\log \frac{|{\Gamma }_{\mathrm{pr}}|}{|{\Gamma }_{\mathrm{pos}}|}.$$

$${\Gamma }_{\mathrm{pos}}={\Gamma }_{\mathrm{pr}}-K{K}^{\top },$$

$$K{K}^{\top }={\Gamma }_{\mathrm{pr}}{G}^{\top }{\Gamma }_{\mathrm{y}}^{-1}G{\Gamma }_{\mathrm{pr}}$$

$${\mathcal{M}}_r=\left\{{\Gamma }_{\mathrm{pr}}-B{B}^{\top }\succ 0:\mathrm{rank}(B)\le r\right\}.$$

$$d_{\mathcal{F}}^2(A,B)=\mathrm{tr}\left[{\ln }^2(A^{-1/2}B{A}^{-1/2})\right]=\sum _{i=1}^n{\ln }^2({\sigma }_i),$$

$$\mathcal{D}(A,B)=\sum _{i=1}^{n}f({\sigma }_i),$$

$${\hat{\Gamma }}_{\mathrm{pos}}=\arg \underset{\Gamma \in \mathcal{M}}{\min }\mathcal{D}({\Gamma }_{\mathrm{pos}},\Gamma ).$$

$$\hat{H}=S_{\mathrm{pr}}^{\top }H{S}_{\mathrm{pr}},$$

$${\hat{\Gamma }}_{\mathrm{pos}}={\Gamma }_{\mathrm{pr}}-B{B}^{\top },B{B}^{\top }=\sum _{i=1}^r{\delta }_i^2{\left(1+{\delta }_i^2\right)}^{-1}{\hat{\mathbf{w}}}_i{\hat{\mathbf{w}}}_i^{\top },{\hat{\mathbf{w}}}_i=S_{\mathrm{pr}}{\mathbf{w}}_i.$$

$$\mathcal{D}({\Gamma }_{\mathrm{pos}},{\hat{\Gamma }}_{\mathrm{pos}})=f(1)r+f\sum _{i=r+1}^{n}f(1/(1+{\delta }_i^2)).$$

$$\hat{L}(\mathbf{\theta },\mathbf{z}):=-\frac{1}{2}{\mathbf{z}}^{\top }{\hat{\Gamma }}_{\mathrm{pos}}\mathbf{z}+\frac{1}{2}\log \frac{|{\Gamma }_{\mathrm{pr}}|}{|{\hat{\Gamma }}_{\mathrm{pos}}|},$$

$${\mathcal{M}}_r^{\prime }=\left\{{\hat{\Gamma }}_{\mathrm{pos}}=({\Gamma }_{\mathrm{pr}}-B{B}^{\top }):{\hat{\Gamma }}_{\mathrm{pos}}-{\Gamma }_{\mathrm{pos}}⪰0,\mathrm{rank}(B)\le r\right\}$$

$$\underset{{\hat{\Gamma }}_{\mathrm{pos}}\in {\mathcal{M}}_r^{\prime }}{\min }\left|\log \frac{|{\Gamma }_{\mathrm{pos}}|}{|{\hat{\Gamma }}_{\mathrm{pos}}|}\right|,$$

$$\log \frac{|{\Gamma }_{\mathrm{pr}}|}{|{\hat{\Gamma }}_{\mathrm{pos}}|}=\sum _{i=1}^r\log (1+{\delta }_i^2)and\log \frac{|{\hat{\Gamma }}_{\mathrm{pos}}|}{|{\Gamma }_{\mathrm{pos}}|}=\sum _{i=r+1}^n\log (1+{\delta }_i^2).$$

$$f(x)=|\log x|.$$

$$\log \frac{|{\Gamma }_{\mathrm{pr}}|}{|{\Gamma }_{\mathrm{pos}}|}=\sum _{i=1}^n\log (1+{\delta }_i^2)={\mathcal{D}}_1({\Gamma }_{\mathrm{pr}},{\Gamma }_{\mathrm{pos}}),$$

$${\mathcal{D}}_1({\Gamma }_{\mathrm{pr}},{\Gamma }_{\mathrm{pos}})={\mathcal{D}}_1({\Gamma }_{\mathrm{pr}},{\hat{\Gamma }}_{\mathrm{pos}})+{\mathcal{D}}_1({\hat{\Gamma }}_{\mathrm{pos}},{\Gamma }_{\mathrm{pos}}),$$

$$\underset{{\hat{\Gamma }}_{\mathrm{pos}}\in {\mathcal{M}}_r^{\prime }}{\min }{\mathcal{D}}_1({\Gamma }_{\mathrm{pos}},{\hat{\Gamma }}_{\mathrm{pos}}).$$

$$\underset{{\hat{\Gamma }}_{\mathrm{pos}}\in {\mathcal{M}}_r^{\prime }}{\min }\underset{\mathbf{z}\in Z_c}{\max }|{\mathbf{z}}^{\top }{\Gamma }_{\mathrm{pos}}\mathbf{z}-{\mathbf{z}}^{\top }{\hat{\Gamma }}_{\mathrm{pos}}\mathbf{z}|$$

$${\mathbf{z}}^{\top }{\hat{\Gamma }}_{\mathrm{pos}}\mathbf{z}={\mathbf{z}}^{\top }{\Gamma }_{\mathrm{pr}}\mathbf{z}-{\mathbf{z}}^{\top }\mathbf{b},\mathbf{b}=\sum _{i=1}^r\frac{{\delta }_i^2}{1+{\delta }_i^2}({\hat{\mathbf{w}}}_i^{\top }\mathbf{z}){\hat{\mathbf{w}}}_i,$$

$$|{\mathbf{z}}^{\top }({\hat{\Gamma }}_{\mathrm{pos}}-{\Gamma }_{\mathrm{pos}})\mathbf{z}|={\mathbf{z}}^{\top }\left(\sum _{i=r+1}^n\frac{{\delta }_i^2}{1+{\delta }_i^2}({\hat{\mathbf{w}}}_i^{\top }\mathbf{z}){\hat{\mathbf{w}}}_i\right).$$

Full text

• March 2018

1 Introduction

Bayesian inference approaches have become increasingly popular as a tool to solve inverse problems [33, 17, 32]. In this setting, the parameters of interest are treated as random variables, endowed with a prior distribution that encodes information available before the data are observed. Observations are modeled by their joint probability distribution conditioned on the parameters of interest, which defines the likelihood function; this distribution incorporates the forward model and a stochastic description of measurement or model errors. The posterior distribution of the unknown parameters is then obtained via Bayes’ formula, which combines the information provided by the prior and the data. Here we focus on a particular class of Bayesian inverse problems, where the forward model is linear and both the prior and the error distributions are Gaussian. Such linear–Gaussian inverse problems arise in many real-world applications [33, 17].

In practice, a difficulty in implementing Bayesian inference for such problems is that the prior and/or the likelihood function may contain unspecified hyperparameters. For instance, the prior distribution usually encodes some kind of correlation among the inversion parameters, but the correlation length is often unknown. Alternatively, the variance and/or the mean of the observational errors may not be available, thus yielding unspecified hyperparameters in the likelihood function. A natural idea to address such problems is to construct a hierarchical Bayesian model for both the inversion parameters and the hyperparameters. Then both sets of parameters can be inferred from data through the characterization of their joint posterior distribution. This inferential approach is often called “fully Bayesian.” A fully Bayesian method can be prohibitively expensive for large-scale inverse problems, however, as it requires sampling from a posterior distribution that is higher dimensional and no longer Gaussian. An alternative strategy is the empirical Bayes (EB) approach [7, 6], which first estimates the hyperparameters by maximizing their marginal likelihood function, i.e., the probability density of the data conditioned only on the hyperparameters, and then plugs in the estimated values to compute the posterior of the inversion parameters. The use of the EB method has been theoretically justified [23, 28]: roughly speaking, these theoretical results show that, under certain conditions and for a sufficiently large sample size, the EB method leads to similar inferential conclusions as fully Bayesian inference. At the same time, the EB method is more computationally efficient than the fully Bayesian approach, as it does not require characterizing the joint posterior distribution.

Nonetheless, the numerical implementation of EB methods remains computationally taxing when dimensionality of the inversion parameters is high, as it involves maximizing the marginal likelihood function, which in turn requires “integrating out” the inversion parameters for each evaluation. (Details about the computational cost of evaluating the marginal likelihood function are given in Section 2.) The goal of this work is to present an approximate EB method that can significantly reduce the computational cost of solving linear inverse problems with hyperparameters. Specifically, our method evaluates the marginal likelihood function by using a low-rank update approximation of the posterior covariance matrix, introduced in previous work for the non-hierarchical setting, e.g., [11, 3, 21, 1, 30, 29]. The intuition behind the low-rank update approximation is that the data may be informative, relative to the prior, only on a low-dimensional subspace of the entire parameter space. As a result, one can consider approximations of the posterior covariance matrix in the form of low-rank negative semidefinite updates of the prior covariance matrix, where the update is obtained by solving a generalized eigenproblem involving the log-likelihood Hessian and the prior precision matrix. The optimality properties of this method are studied in [30], which shows that this low-rank approximation is optimal with respect to a specific class of loss functions.

Using the approximation of the posterior covariance developed in [11, 3, 30], we introduce a new method to efficiently compute the marginal likelihood function. We prove that the proposed method yields an optimal approximation of the marginal likelihood function in a minimax sense. Our last contribution lies in the numerical implementation. Unlike the inverse problems considered in the work mentioned above, where the prior is fixed, our problems require repeatedly computing the square roots of different prior covariance matrices—which can be prohibitively expensive when the dimensionality of the problem is high. To address the issue, we use the spectral approximation developed in [15] to compute the square root of the prior covariance matrix, resulting in a very efficient algorithm for evaluating (and maximizing) the marginal likelihood function.

The rest of this paper is organized as follows. In Section 2 we introduce the basic setup of the empirical Bayesian method for solving inverse problems. The low-rank update approximation for evaluating the log marginal likelihood function is developed and analyzed in Section 3, and the detailed numerical implementation of the proposed method is discussed in Section 4. Section 5 provides two numerical examples to demonstrate the performance of the proposed method, and Section 6 offers concluding remarks.

2 Problem setup

Consider a linear inverse problem in a standard setting:

[TABLE]

where $\mathbf{y}\in{\mathbb{R}}^{m}$ is the data, $\mathbf{x}\in{\mathbb{R}}^{n}$ is the unknown, $G$ is a $n\times m$ matrix, often called the forward operator, and $\bm{\eta}$ is the observation noise. Such an inverse problem can be solved with a Bayesian inference method: namely we assume that the prior probability density of $\mathbf{x}$ is $\pi(\mathbf{x})$ and the likelihood function is $\pi(\mathbf{y}|\mathbf{x})$, and thus the posterior is given by Bayes’ theorem:

[TABLE]

Throughout this paper, we use the notation $\pi(\cdot)$ as a generic representation of probability density; the specific density is made clear by its arguments. We further assume a Gaussian likelihood and a Gaussian prior with a non-singular covariance matrix $\Gamma_{\rm pr}\succ 0$ and, without loss of generality, zero mean:

[TABLE]

In this setting, the resulting posterior distribution is also Gaussian

[TABLE]

with mean and covariance matrix given by

[TABLE]

where

[TABLE]

is the Hessian of the log-likelihood.

Now we consider a more complex situation, where the matrices $G$, $\Gamma_{\rm pr}$ and $\Gamma_{\rm obs}$ (or some of them) depend on a vector of unspecified hyperparameters $\bm{\theta}\in{\mathbb{R}}^{p}$. As mentioned in the introduction, a popular method to deal with such problems is the empirical Bayesian (EB) approach, which determines $\bm{\theta}$ by maximizing the marginal likelihood function:

[TABLE]

or the marginal posterior density:

[TABLE]

where $p(\bm{\theta})$ is the prior density of $\bm{\theta}$. Note that the computations of these two objectives are very similar and so we shall only discuss (4). It is easy to see that the optimization problem in (4) is equivalent to minimizing the negative log marginal likelihood:

[TABLE]

where $G$, $\Gamma_{\rm pr}$ and $\Gamma_{\rm obs}$ depend on $\bm{\theta}$ and $\mathbf{z}=G^{T}\Gamma_{\rm obs}^{-1}\mathbf{y}$. Note that an important special case arises when only the prior depends on the hyperparameter $\bm{\theta}$ and the likelihood function is fixed; in this case, we can simply minimize

[TABLE]

For conciseness, below we will present our method for (6), while noting that all the results can be trivially extended to (5).

Direct evaluation of (6) is not desirable for large scale problems, as it requires several operations with $O(n^{3})$ complexity. In what follows, we present an accurate and efficient—with $O(n^{2}r)$ complexity for some $r\ll n$—method to approximate $L(\bm{\theta})$, based on a rank-$r$ update approximation of $\Gamma_{\rm pos}$.

3 The optimal approximation method

3.1 Optimal low-rank update approximation of the posterior covariance

The proposed method begins with the optimal low-rank update approximation of posterior covariance developed in [3, 30], which is briefly reviewed here. Note that $\Gamma_{\rm pos}$ can be written as a non-positive update of $\Gamma_{\rm pr}$:

[TABLE]

where

[TABLE]

and $\Gamma_{\rm y}=\Gamma_{\rm obs}+G\,\Gamma_{\rm pr}\,G^{\top}$ is the covariance of the marginal distribution of $\mathbf{y}$. This update of $\Gamma_{\rm pr}$ is negative semidefinite because the data add information; they cannot increase the prior variance in any direction. As has been discussed in previous work [30], in many practical problems, the low-rank structure often lies in the update of $\Gamma_{\rm pr}$ that yields $\Gamma_{\rm pos}$, rather than in $\Gamma_{\rm pos}$ itself. Hence, [30] proposes to use the following class of positive definite matrices to approximate $\Gamma_{\rm pos}$:

[TABLE]

To establish optimality statements regarding the approximation of a covariance matrix, [30] adopts the following distance between symmetric positive definite (SPD) matrices of size $n$:

[TABLE]

where $(\sigma_{i})_{i=1}^{n}$ are the generalized eigenvalues of the pencil $(A,B)$. This metric was first introduced by Rao in [26], and can measure the difference between two Gaussian distributions with the same mean. We direct interested readers to [29] for a detailed discussion and other applications of this metric. It is important to note that this distance is generalized in [30] to be

[TABLE]

where $f$ is a function in $C^{1}({\mathbb{R}}^{+})$ that satisfies $f^{\prime}(x)(1-x)<0$ for any $x\neq 1$, and $\lim_{x\rightarrow\infty}f(x)=\infty$. This generalization will be used for the proof of our optimality statement in next section. Thus we seek a low-rank approximation $\widehat{\Gamma}_{\rm pos}$ of the covariance matrix $\Gamma_{\rm pos}$, such that

[TABLE]

The solution of (8) can be derived analytically regardless of the specific choice of $f$. Specifically, let $S_{\rm pr}$ be any square root of the prior covariance matrix such that $\Gamma_{\rm pr}=S_{\rm pr}\,S_{\rm pr}^{\top}$. We define the prior-preconditioned Hessian as

[TABLE]

which plays an essential role in our numerical method. Now let $(\delta_{i}^{2},\mathbf{w}_{i})$ be the eigenvalue-eigenvector pairs of $\widehat{H}$ with the ordering $\delta_{i}^{2}\geq\delta_{i+1}^{2}$; then a solution of (8) is given by:

[TABLE]

The corresponding minimum distance is

[TABLE]

The minimizer (10) is unique if the first $r$ eigenvalues of $\widehat{H}$ are distinct.

3.2 Approximating the log-likelihood function

Now we apply the optimal low rank approximation to our problem. The idea is rather straightforward: we approximate $L$ in (6) with

[TABLE]

for some approximate posterior covariance matrix $\widehat{\Gamma}_{\rm pos}$. Next we shall derive the appropriate choice of $\widehat{\Gamma}_{\rm pos}$. To do this, we need to impose an additional assumption on the approximate posterior covariance matrix: $\widehat{\Gamma}_{\rm pos}-\Gamma_{\rm pos}\succeq 0$, which means that the approximation itself should not create any new information. As a result, the class of matrices for approximating $\Gamma_{\rm pos}$ becomes

[TABLE]

for some maximum rank $r$. Next we shall consider the approximations of the two terms $\log({|\Gamma_{\rm pr}|}/{|\Gamma_{\rm pos}|})$ and $\mathbf{z}^{\top}\,\Gamma_{\rm pos}\,\mathbf{z}$ in (6) separately.

3.2.1 Approximating the log–determinant term.

First, we consider finding a matrix $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$ to approximate $\log({|\Gamma_{\rm pr}|}/{|\Gamma_{\rm pos}|})$ with $\log({|\Gamma_{\rm pr}|}/{|\widehat{\Gamma}_{\rm pos}|})$. In this setting, it is easy to see that the approximation error in the log-marginal likelihood (6) is $|\log({|\Gamma_{\rm pos}|}/{|\widehat{\Gamma}_{\rm pos}|})|$, and a natural way to determine $\widehat{\Gamma}_{\rm pos}$ is to seek a $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$ that minimizes this approximation error. To this end, we have the following theorem:

Theorem 3.1

Suppose that we approximate $\log({|\Gamma_{\rm pr}|}/{|\Gamma_{\rm pos}|})$ with $\log({|\Gamma_{\rm pr}|}/{|\widehat{\Gamma}_{\rm pos}|})$ for some $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$. The matrix $\widehat{\Gamma}_{\rm pos}\in{\cal M}^{\prime}_{r}$ that minimizes the resulting approximation error, i.e., the solution of

[TABLE]

is given by (10). Moreover, the optimal approximation and the associated error are, respectively,

[TABLE]

• Proof

We prove this theorem using the optimality results in [30]. To start, we choose a particular distance metric by letting

[TABLE]

We denote the resulting distance metric as ${\mathcal{D}}_{1}(A,B)$ to indicate that such a metric is actually the $1$-norm of $(\ln(\sigma_{1}),\,\ldots,\,\ln(\sigma_{n}))$ while $d_{{\cal F}}$ is the $2$-norm of the same sequence. It can be verified that

[TABLE]

as $\log(1+\delta_{i}^{2})\geq 0$ for all $i\in{\mathbb{N}}$. Since the approximate posterior covariance $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$, we can show

[TABLE]

where ${\mathcal{D}}_{1}(\Gamma_{\rm pr},\widehat{\Gamma}_{\rm pos})=\log({|\Gamma_{\rm pr}|}/{|\widehat{\Gamma}_{\rm pos}|})$ is the approximation and ${\mathcal{D}}_{1}(\widehat{\Gamma}_{\rm pos},\Gamma_{\rm pos})=\log({|\widehat{\Gamma}_{\rm pos}|}/{|\Gamma_{\rm pos}|})$ is the error associated with it. Note that (15) does not hold without the assumption $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$. Thus (13) can be rewritten as,

[TABLE]

Recall that the solution of $\min_{\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}}{\mathcal{D}}_{1}(\Gamma_{\rm pos},\widehat{\Gamma}_{\rm pos})$ is given by (10), and it is easy to verify that the matrix $\widehat{\Gamma}_{\rm pos}$ given by (10) is in ${\cal M}_{r}^{\prime}$, which implies that (10) also provides the solution of (16). As a result, the optimal approximation and the associated approximation error are given by (14), which completes the proof.

3.2.2 Approximating the quadratic term.

Similarly, we can also find an approximate posterior covariance $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$ and approximate $\mathbf{z}^{\top}\,\Gamma_{\rm pos}\,\mathbf{z}$ with $\mathbf{z}^{\top}\,\widehat{\Gamma}_{\rm pos}\,\mathbf{z}$. This problem is a bit more complicated: since $\mathbf{z}$ (which is a linear transformation of the data $\mathbf{y}$) is random, we cannot determine the matrix $\widehat{\Gamma}_{\rm pos}$ by directly minimizing the approximation error. In this case, a useful treatment is to apply the minimax criterion, i.e., to seek a matrix $\widehat{\Gamma}_{\rm pos}\in{\cal M}^{\prime}_{r}$ that minimizes the maximum approximation error with respect to $\mathbf{z}$. In particular, for the maximum error to exist, we shall require $\mathbf{z}$ to be bounded: $\mathbf{z}\in Z_{c}=\{\mathbf{z}:\|\mathbf{z}\|_{2}\leq c\}$ for a constant $c>0$. (See Remark 3.1 for a discussion of this boundedness assumption.) The following theorem provides the optimal solution to this problem.

Theorem 3.2

Suppose that we approximate $\mathbf{z}^{\top}\,\Gamma_{\rm pos}\,\mathbf{z}$ with $\mathbf{z}^{\top}\,\widehat{\Gamma}_{\rm pos}\,\mathbf{z}$, for some $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$. The matrix $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$ that achieves the minimax approximation error, i.e., the solution of

[TABLE]

is given by (10). Moreover, the resulting approximation is

[TABLE]

*and the associated approximation error is *

[TABLE]

• Proof

For any given $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$, it is easy to see that the solution of

[TABLE]

is $\mathbf{z}=c{\mathbf{v}}_{\max}$ where $\mathbf{v}_{\max}$ is the eigenvector of the largest eigenvalue of $\widehat{\Gamma}_{\rm pos}-\Gamma_{\rm pos}$, and the maximum error is $c^{2}\rho(\widehat{\Gamma}_{\rm pos}-\Gamma_{\rm pos})$ where $\rho(\cdot)$ is the spectral radius. Thus (17) becomes $\min_{\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}}\rho(\widehat{\Gamma}_{\rm pos}-\Gamma_{\rm pos})$, or equivalently,

[TABLE]

and it follows immediately that the optimal $\widehat{\Gamma}_{\rm pos}$ is given by (10). Substituting (10) into $\mathbf{z}^{\top}\,\widehat{\Gamma}_{\rm pos}\,\mathbf{z}$ and $\mathbf{z}^{\top}\,(\widehat{\Gamma}_{\rm pos}-\Gamma_{\rm pos})\,\mathbf{z}$ yields (18) and (19), respectively.

In principle, of course, $\|\bm{z}\|_{2}$ is not bounded from above, since $\bm{z}$ is normally distributed. However, imposing the boundedness assumption, which considerably simplifies the theoretical analysis, does not limit the applicability of the method or affect the optimal solution. Specifically, we have the following remarks:

Remark 3.1

First, since $\bm{z}$ follows a Gaussian distribution, one can always choose a constant $c$ such that the inequality $\|\bm{z}\|_{2}\leq c$ holds with probability arbitrarily close to one. Second, the minimax solution $\widehat{\Gamma}_{\rm pos}$ is independent of the value of $c$.

Now we combine the two approximate treatments, which is essentially to approximate (6) by (11) with $\widehat{\Gamma}_{\rm pos}$ given by (10). It is easy to see that the approximation error is

[TABLE]

and we have the following result regarding its optimality:

Corollary 3.3

Suppose that we approximate (6) with (11) for some matrix $\widehat{\Gamma}_{\rm pos}\in{\cal M}_{r}^{\prime}$. The matrix $\widehat{\Gamma}_{\rm pos}$ given by (10) achieves the minimax approximation error, i.e., it solves

[TABLE]

As both terms on the right hand side of (20) are nonnegative, the corollary is a direct consequence of Theorems 3.1 and 3.2.

4 Numerical implementation

Here we discuss the numerical implementation of our approximate method to evaluate (6). In principle, this involves two computationally intensive components—both requiring $O(n^{3})$ computations under standard numerical treatments. The first task is to compute the eigenvalues of $\widehat{H}$. Recall that our method only requires the $r$ leading eigenvalues and associated eigenvectors of $\widehat{H}$, which can be computed efficiently with a randomized algorithm for the singular value decomposition (SVD) [20, 13]. The second task is to compute $S_{\rm pr}$, the square root of $\Gamma_{\rm pr}$. As will become clear later, we do not necessarily need $S_{\rm pr}$; rather, our algorithm only requires the ability to evaluate $S_{\rm pr}\Omega$ for a given matrix $\Omega$. To this end, we resort to the approximation method proposed in [15]. We provide a brief description of the two adopted methods here, tailored according to our specific purposes.

4.1 Randomized SVD

The main idea behind the randomized SVD approach is to identify a subspace that captures most of the action of a matrix, using random sampling [20, 13, 10]. The original matrix is then projected onto this subspace, yielding a smaller matrix, and a standard SVD is then applied to the smaller matrix to obtain the leading eigenvalues of the original matrix. A randomized algorithm for constructing the rank–$r^{\prime}$ SVD of an $n\times n$ matrix $\widehat{H}$, given in [13], proceeds as follows:

Draw an $n\times r^{\prime}$ Gaussian random matrix $\Omega$; 2. 2.

Form the $n\times r^{\prime}$ sample matrix $Y=\widehat{H}\Omega$; 3. 3.

Form an $n\times r^{\prime}$ matrix $Q$ with orthonormal columns, such that $Y=QR$; 4. 4.

Form the $r^{\prime}\times n$ matrix $B=Q^{T}\widehat{H}$; 5. 5.

Compute the SVD of the $r^{\prime}\times n$ matrix $B=\hat{U}\Sigma V^{T}$.

According to [13], for this algorithm to be robust, it is important to oversample a bit, namely to choose $r^{\prime}>r$ if $r$ is the desired rank. One can obtain a probabilistic error bound, i.e., that

[TABLE]

holds with a probability at least $1-6(r^{\prime}-r)^{-(r^{\prime}-r)}$, under some very mild assumptions on $r^{\prime}$ [13]. Note that the diagonal entries of $\Sigma$ are the eigenvalues of $\widehat{H}$ and the columns of $V$ are the eigenvectors. Since we obtain $r^{\prime}$ eigenpairs in the algorithm, we take $\{\delta_{i},\mathbf{w}_{i}\}_{i=1}^{r}$ to be the $r$ dominant ones among them. Above we only present the basic implementation of the randomized SVD method; further details and possible improvements of the method can be found in [13] and the references therein. Finally we also want to emphasize that our approximate EB method can use the eigenvalues/eigenvectors computed with any numerical approach; it is not tailored to or tied with the randomized SVD.

4.2 The Chebyshev spectral method for computing $S_{\rm pr}$

Recall that $\widehat{H}=S_{\rm pr}^{\top}H\,S_{\rm pr}$, and as a result the randomized SVD only requires evaluating ${S_{\rm pr}}\Omega$ where $\Omega$ is a randomly generated $n\times r^{\prime}$ matrix. Here we adopt the approximate method proposed in [15], which is based on the following lemma:

Lemma 4.1** (Lemma 2.1 in [15])**

Suppose that $D$ is a real symmetric positive definite matrix. Then there exists a polynomial $p(\cdot)$ such that $\sqrt{D}=p(D)$, and the degree of $p$ is equal to the number of distinct eigenvalues of $D$ minus $1$.

Though $\sqrt{D}$ is exactly equal to a polynomial of $D$, the degree of the polynomial might be very large if $D$ has a large number of distinct eigenvalues. Thus, instead of trying to find the exact polynomial $p$ in $D$ that equals $\sqrt{D}$, the aforementioned work computes a Chebyshev approximation to it.

It is clear that $D=O\Lambda O^{T}$, where $\Lambda=\mathrm{diag}(\lambda_{1}\cdots\lambda_{n})$, and $\lambda_{1}\geq\cdots\geq\lambda_{n}$. Now suppose that we have a polynomial $p(\cdot)$ such that $p(\lambda)\approx\sqrt{\lambda}$. Then we have

[TABLE]

That is, once we have an approximation of ${\sqrt{\lambda}}$ in the interval $[\lambda_{\min},\lambda_{\max}]$, where $\lambda_{\min}$ and $\lambda_{\max}$ are respectively lower and upper bounds on the eigenvalues of $D$, we immediately get an approximation of $\sqrt{D}$. A popular method to construct the approximation $p(\cdot)$ is the Chebyshev polynomial interpolation. The standard Chebyshev polynomials are defined in $[-1,\,1]$ as [24],

[TABLE]

and the associated Chebyshev points are given by

[TABLE]

It is well known that the Chebyshev interpolant has spectral accuracy for analytical functions on $[-1,1]$; see, e.g. [12].

Here since we intend to approximate the function $\sqrt{\lambda}$ over the range $[\lambda_{\min},\lambda_{\min}]$ rather than $[-1,\,1]$, we shall use the scaled and shifted Chebyshev polynomials:

[TABLE]

and the associated Chebyshev points become

[TABLE]

The interpolant $p_{k}(x)$ can be expressed as [24],

[TABLE]

where the coefficients $c_{i}$ are given by

[TABLE]

It is easy to verify that the scaled and shifted Chebyshev polynomials satisfy the following recurrence relation,

[TABLE]

with

[TABLE]

Now recall that we actually want to compute $\sqrt{D}\Omega$. Taking advantage of the recurrence relation (27), we obtain,

[TABLE]

where

[TABLE]

for $i=1,\ldots,k-1$. The complete procedure for computing $\sqrt{D}\Omega$ with the Chebyshev approximation is given in Algorithm 1.

Finally we provide some remarks regarding the implementation of the method:

• •

The original algorithm present in [15] is to compute the product of $\sqrt{D}$ and a vector, but as is shown in Algorithm 1, its extension to the computation of the product of $\sqrt{D}$ and a matrix is straightforward.

• •

The method requires upper and lower bounds on the eigenvalues of $D$, which are computed with an algorithm based on the safeguarded Lanczos method [15]. As is discussed in [15], these estimates need not be of high accuracy.

• •

The error bound of the proposed Chebyshev approximation is given by Theorems 3.3 and 3.4 in [15].

• •

If desired, all the matrix-vector multiplications can be performed with the fast multipole method to further improve efficiency [15].

4.3 The complete algorithm

Now we summarize the complete scheme for constructing the approximation $\widehat{L}(\bm{\theta},\mathbf{z})$ of (6) with a given rank $r$, using the methods discussed in Sections 4.1 and 4.2:

Compute the first $r$ eigenpairs, $\{(\delta_{i}^{2},\mathbf{w}_{i})\}_{i=1}^{r}$, of $\widehat{H}(\bm{\theta})=S_{\rm pr}HS_{\rm pr}$, using the randomized SVD method and the Chebyshev spectral method. 2. 2.

Let $\widehat{\mathbf{w}}_{i}=S_{\rm pr}\mathbf{w}_{i}$ and evaluate

[TABLE]

5 Examples

5.1 An image deblurring problem

We first test our method on an imaging deblurring problem, which involves recovering a latent image from noisy observations of a blurred version of the image [14]. In particular, it is assumed that the blurred image is obtained as a convolution of the latent image with a point spread function (PSF), and as a result the forward model is:

[TABLE]

where $f_{\mathrm{PSF}}(t_{1},t_{2})$ is the PSF and $D$ is the domain of the image. In this example we take the image domain to be $D=[-1,1]^{2}$ and the PSF to be

[TABLE]

where $t$ is a parameter controlling the size of the spreading. Moreover, we assume that the data $y$ are measured at $m=64^{2}=4096$ observation locations evenly distributed in $D$, and that the observation errors are mutually independent and Gaussian with zero mean and variance $(0.1)^{2}$. We represent the unknown $x$ on $256\times 256$ mesh grid points, and thus the dimensionality of the inverse problem is $n=65536$. The prior on $x$ is a Gaussian distribution with zero mean and covariance kernel [31, 27]:

[TABLE]

where $d=\|\mathbf{t}-\mathbf{t}^{\prime}\|_{2}$, $\Gamma(\cdot)$ is the Gamma function, and $B_{\nu}(\cdot)$ is the modified Bessel function. (30) is known as the Matérn covariance, and several authors have suggested that such covariances can often provide better models for many realistic physical processes than the popular squared exponential covariance [31, 27]. A random function with the Matérn covariance is $[\nu-1]$ mean-square (MS) differentiable. In this example, we choose $\nu=3$ implying second order MS differentiability. We set the standard deviation $\sigma$ to one. The correlation length $\rho$ is treated as a hyperparameter to be inferred in this example.

We now use the proposed empirical Bayesian method to solve the inverse problem. We first assume that the true correlation length is $\rho=1$, and randomly generate a true image from the associated prior, shown in Fig. 1 (left). We then test two cases of the forward problem, with $t=0.002$ and $t=0.02$, where the latter yields an inverse problem that is much more ill-posed than the former. We assume that the data are observed on a $128\times 128$ uniformly distributed mesh. We then apply the two convolution operators to the generated image and add observational noise to the results, producing the synthetic data also shown in Fig. 1.

Using the proposed approximate EB method, we evaluate the negative log marginal likelihood function $L$ over a range of $\rho$ values, for both cases of $t$, and plot the results in Fig 2. For $t=0.002$ (left figure), we compute $L$ with ranks $r=300$, $400$ and $500$, and $600$, and observe that the curves converge as the rank increases; indeed, the results with $r=500$ and $600$ appear identical. For $t=0.02$—because the problem is more ill-posed—we can implement the method with ranks $r=50$, $75$, $100$, and $150$ and observe convergence. In particular, while the results of $r=50$ deviate from the others, the results with $r=75$, $100$, and $150$ look nearly identical, implying that $r=75$ is sufficient for an accurate approximation of the marginal likelihood in this case. For both cases, the optimal value of $\rho$ is found to be $0.1$, which is actually the true hyperparameter value. We then compute the posterior mean of $x$, after fixing $\rho=0.1$ in the Gaussian prior on $x$, and show the results in Fig. 3.

Note that the EB method is able to find the appropriate hyperparameter values in this example. Nevertheless, our intention is not to illustrate that the EB method can always identify the correct value of the hyperparameters; rather, the main purpose of the example is to demonstrate that, should one choose to use the EB method, the proposed approximations can efficiently and accurately compute the marginal likelihood.

5.2 X-ray computed tomography

Our second example is an X-ray computed tomography (CT) problem. X-ray CT is a popular diagnostic tool for medical tomographic imaging of the human body. It provides images of the object by assigning an X-ray attenuation coefficient to each pixel [22]. Specifically let $u$ denote the unknown image supported in the unit ball $B({\mathbf{0}},1)$ in ${\mathbb{R}}^{2}$. In the case of two-dimensional parallel beam CT, the projection data (or sinogram) $f$ for each $\varphi\in[0,2\pi)$ and $s\in{\mathbb{R}}$ is obtained via the following Radon transform [25]:

[TABLE]

where $\bm{\theta}=(\cos\varphi,\sin\varphi)$ and $\bm{\theta}^{\perp}=(-\sin\varphi,\cos\varphi)$. We test the problem using a $256\times 256$ ground truth image shown in Fig. 5 (left), which is taken from the Harvard Whole Brain Atlas [16]. The projection data, shown in Fig. 5 (right), is simulated by plugging the true image into the Radon transform and adding measurement noise to the result. The noise is taken to be Gaussian with zero mean and a signal-to-noise ratio of $1\%$. The size of the discrete Radon transform operator depends on both the size of the image and the number of projections. In our simulation, we used $60$ projections equispatially sampled from [math] to $\pi$. To avoid an inverse crime [18], we recover the image with a $128\times 128$ resolution, and thus the dimensionality of the inverse problem is $128^{2}$. We take the prior distribution to be zero-mean Gaussian, with a covariance matrix given by an anisotropic Matérn kernel:

[TABLE]

where

[TABLE]

Now we consider all four parameters $\nu$, $\sigma$, $\rho_{1}$ and $\rho_{2}$ as hyperparameters to be determined, and we also assume that the variance $\epsilon^{2}$ of the measurement noise is not available. We will identify the five parameters using the proposed low-rank approximate EB method.

We solve the optimization problems with 9 different ranks ($r=500$, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500), and plot the optimal values of the five hyperparameters and the associated value of $L$ versus the rank $r$ in Fig. 6. As one can see, the obtained optimal values of the hyperparameters converge as the rank increases, and in particular $r=4000$ is sufficient to obtain accurate estimates of the hyperparameters. Next we shall illustrate that the obtained hyperparameter values can lead to rather good inference results. To show this, we compute the posterior mean and the variance of the image using the hyperparameter values computed with the various ranks, and we then compute the peak noise-to-signal ratio (PSNR) of the posterior means against the true image, which is a commonly used measure of quality for image reconstruction. Note that the low rank approximation is only used to obtain the hyperparameters; given these hyperparameter values, the posterior means and variances are computed directly without approximation. We plot the PSNR as a function of the rank $r$ in Fig. 7. We also show the posterior mean and variance field associated with each rank in the plot. One can see from the figure that, as the rank increases, the resulting hyperparameters can recover images with higher quality in terms of PSNR. Finally, in Fig. 8, we plot the posterior mean and variance using the hyperparameters computed with $r=4500$, which are taken to be the inference results of this problem.

These results demonstrate that the outcome of Bayesian inference can depend critically on the hyperparameters, and that, at least in this example, the proposed method can provide reasonable values for these parameters, provided that a sufficiently large rank is used. We also note here that, to better recover images with sharp edges, one might resort to more complex Gaussian hypermodels [2, 4, 5] than the one considered here, or even non-Gaussian models such as the Besov [9, 19] and TV-Gaussian [34] priors. Nevertheless, if one has chosen to use a Gaussian hypermodel, our method can efficiently determine the values of the associated hyperparameters.

6 Conclusions

This paper investigates empirical Bayesian methods for determining hyperparameters in linear inverse problems. We develop an efficient numerical method to approximately evaluate the marginal likelihood function of the hyperparameters for large-scale problems, based on a low-rank approximation of the update from the prior to the posterior covariance. The proposed method can achieve a computational complexity of $O(n^{2}r)$ for $r\ll n$, while a standard full-rank or direct evaluation approach requires computations of $O(n^{3})$ complexity. We also show that this approximation of the marginal likelihood is optimal in a minimax sense. Two numerical examples illustrate that the proposed algorithm can accurately evaluate the marginal likelihood function required in the EB method. This approach may be useful in a wide range of applications where the unknown is of very high dimension, such as medical and geophysical image reconstruction. Finally, it is worth noting that, while the approximate EB method presented in this work can only be applied to linear inverse problems, there have been several efforts successfully extending low-rank posterior approximations to nonlinear problems [21, 8]; along these lines, we expect that our method could also be extended to EB computations for nonlinear inverse problems.

Acknowledgement

Li was partially supported by the NSFC under grant number 11771289. Marzouk was partially supported by the US Department of Energy, Office of Advanced Scientific Computing Research, under grant number DE-SC0009297.