Bayesian inference and mathematical imaging. Part III: probability & convex optimisation. Dr. Marcelo Pereyra http://www.macs.hw.ac.uk/ ∼ mp71/ Maxwell Institute for Mathematical Sciences, Heriot-Watt University January 2019, CIRM, Marseille. M. Pereyra (MI — HWU) Bayesian mathematical imaging 0 / 53
Outline Bayesian inference in imaging inverse problems 1 MAP estimation with Bayesian confidence regions 2 Posterior credible regions Uncertainty visualisation Hypothesis testing A decision-theoretic derivation of MAP estimation 3 4 Hierarchical MAP estimation with unknown regularisation parameters Conclusion 5 M. Pereyra (MI — HWU) Bayesian mathematical imaging 1 / 53
Imaging inverse problems We are interested in an unknown image x ∈ R d . We measure y , related to x by a statistical model p ( y ∣ x ) . The recovery of x from y is ill-posed or ill-conditioned, resulting in significant uncertainty about x . For example, in many imaging problems y = Ax + w , for some operator A that is rank-deficient, and additive noise w . M. Pereyra (MI — HWU) Bayesian mathematical imaging 2 / 53
The Bayesian framework We use priors to reduce uncertainty and deliver accurate results. Given the prior p ( x ) , the posterior distribution of x given y p ( x ∣ y ) = p ( y ∣ x ) p ( x )/ p ( y ) models our knowledge about x after observing y . In this talk we consider that p ( x ∣ y ) is log-concave; i.e., p ( x ∣ y ) = exp {− φ ( x )}/ Z , where φ ( x ) is a convex function and Z = ∫ exp {− φ ( x )} d x . M. Pereyra (MI — HWU) Bayesian mathematical imaging 3 / 53
Maximum-a-posteriori (MAP) estimation The predominant Bayesian approach in imaging is MAP estimation x MAP = argmax p ( x ∣ y ) , ˆ x ∈ R d (1) = argmin φ ( x ) , x ∈ R d computed efficiently, even in very high dimensions, by (proximal) convex optimisation (Green et al., 2015; Chambolle and Pock, 2016). However, MAP estimation has some limitations, e.g., 1 it provides little information about p ( x ∣ y ) , 2 it is not theoretically well understood (yet), 3 it struggles with unknown/partially unknown models. M. Pereyra (MI — HWU) Bayesian mathematical imaging 4 / 53
Illustrative example: astronomical image reconstruction Recover x ∈ R d from low-dimensional degraded observation y = M F x + w , where F is the continuous Fourier transform, M ∈ C m × d is a measurement mask operator, and w is Gaussian noise. We use the model p ( x ∣ y ) ∝ exp (−∥ y − M F x ∥ 2 / 2 σ 2 − θ ∥ Ψ x ∥ 1 ) 1 R n + ( x ) . (2) ˆ x MAP ∣ y ∣ Figure : Radio-interferometric image reconstruction of the W28 supernova . M. Pereyra (MI — HWU) Bayesian mathematical imaging 5 / 53
Outline Bayesian inference in imaging inverse problems 1 MAP estimation with Bayesian confidence regions 2 Posterior credible regions Uncertainty visualisation Hypothesis testing A decision-theoretic derivation of MAP estimation 3 4 Hierarchical MAP estimation with unknown regularisation parameters Conclusion 5 M. Pereyra (MI — HWU) Bayesian mathematical imaging 6 / 53
Outline Bayesian inference in imaging inverse problems 1 MAP estimation with Bayesian confidence regions 2 Posterior credible regions Uncertainty visualisation Hypothesis testing A decision-theoretic derivation of MAP estimation 3 4 Hierarchical MAP estimation with unknown regularisation parameters Conclusion 5 M. Pereyra (MI — HWU) Bayesian mathematical imaging 7 / 53
Posterior credible regions Where does the posterior probability mass of x lie? Recall that C α is a ( 1 − α ) % posterior credible region if P [ x ∈ C α ∣ y ] = 1 − α, and the decision-theoretically optimum is the HPD region (Robert, 2001) C ∗ α = { x ∶ φ ( x ) ≤ γ α } , with γ α ∈ R chosen such that ∫ C ∗ α p ( x ∣ y ) d x = 1 − α holds. We could estimate C ∗ α by MCMC sampling, but in high-dimensional log-concave models this is not necessary because something beautiful happens... M. Pereyra (MI — HWU) Bayesian mathematical imaging 8 / 53
A concentration phenomenon! Figure : Convergence to “typical” set { x ∶ log p ( x ∣ y ) ≈ E [ log p ( x ∣ y )]} . M. Pereyra (MI — HWU) Bayesian mathematical imaging 9 / 53
Proposed approximation of C ∗ α Theorem 2.1 (Pereyra (2016)) Suppose that the posterior p ( x ∣ y ) = exp {− φ ( x )}/ Z is log-concave on R d . Then, for any α ∈ ( 4exp (− n / 3 ) , 1 ) , the HPD region C ∗ α is contained by √ C α = { x ∶ φ ( x ) ≤ φ ( ˆ ˜ x MAP ) + d τ α + d )} , √ with universal positive constant τ α = 16log ( 3 / α ) independent of p ( x ∣ y ) . Remark 1 : ˜ C α is a conservative approximation of C ∗ α , i.e., x ∉ ˜ C α � ⇒ x ∉ C ∗ α . Remark 2 : ˜ C α is available as a by-product in any convex inference problem that is solved by MAP estimation! M. Pereyra (MI — HWU) Bayesian mathematical imaging 10 / 53
Approximation error bounds Is ˜ C α a reliable approximation of C ∗ α ? Theorem 2.2 (Finite-dimensional error bound (Pereyra, 2016)) √ γ α = φ ( ˆ x MAP ) + d τ α + d . If p ( x ∣ y ) is log-concave on R d , then Let ˜ γ α − γ α 0 ≤ ˜ ≤ 1 + η α d − 1 / 2 , d √ √ with universal positive constant η α = 16log ( 3 / α ) + 1 / α . C α is stable (as d becomes large, the error ( ˜ γ α − γ α )/ d ⪅ 1). Remark 3: ˜ Remark 4: The lower and upper bounds are asymptotically tight w.r.t. d . M. Pereyra (MI — HWU) Bayesian mathematical imaging 11 / 53
Outline Bayesian inference in imaging inverse problems 1 MAP estimation with Bayesian confidence regions 2 Posterior credible regions Uncertainty visualisation Hypothesis testing A decision-theoretic derivation of MAP estimation 3 4 Hierarchical MAP estimation with unknown regularisation parameters Conclusion 5 M. Pereyra (MI — HWU) Bayesian mathematical imaging 12 / 53
Example: Uncertainty visualisation in astro-imaging Radio-interferometry with redundant wavelet frame (Cai et al., 2017). dirty image ˆ x penML ( y ) ˆ x MAP approx. credible intervals (scale 10 × 10 pixels) “exact” intervals (MCMC, minutes) M31 radio galaxy (size 256 × 256 pixels, comp. time 1 . 8 secs.) M. Pereyra (MI — HWU) Bayesian mathematical imaging 13 / 53
Example: Uncertainty visualisation in astro-imaging Radio-interferometry with redundant wavelet frame (Cai et al., 2017). dirty image ˆ x penML ( y ) ˆ x MAP approx. credible intervals (scale 10 × 10 pixels) “exact” intervals (MCMC, minutes) 3C2888 radio galaxy (size 256 × 256 pixels, comp. time 1 . 8 secs.) M. Pereyra (MI — HWU) Bayesian mathematical imaging 14 / 53
Example: uncertainty visualisation in microscopy Live cell microscopy sparse super-resolution (Zhu et al., 2012): y = Ax + w x MAP = argmin x ∈ R d ∥ y − Hx ∥ 2 / 2 σ 2 + λ ∥ x ∥ 1 ˆ ( A is a blur operator) (log-scale) Consider the molecular structure in the highlighted region: Are we confident about this structure (its presence, position, etc.)? Idea : use ˜ C α to explore/quantify the uncertainty about this structure. M. Pereyra (MI — HWU) Bayesian mathematical imaging 15 / 53
Example: uncertainty visualisation in microscopy Position uncertainty quantification Find maximum molecule displacement within ˜ C 0 . 05 : Molecule position uncertainty x MAP = argmin x ∈ R d ∥ y − Ax ∥ 2 / 2 σ 2 + λ ∥ x ∥ 1 ˆ ( ± 93 nm × ± 140 nm ) Note: Uncertainty analysis ( ± 93 nm × ± 140 nm ) in agreement with MCMC estimates ( ± 78 nm × ± 125 nm - approx. error of order of 1 pixel), and with experimental results (average precision ± 80 nm ) of Zhu et al. (2012). M. Pereyra (MI — HWU) Bayesian mathematical imaging 16 / 53
Outline Bayesian inference in imaging inverse problems 1 MAP estimation with Bayesian confidence regions 2 Posterior credible regions Uncertainty visualisation Hypothesis testing A decision-theoretic derivation of MAP estimation 3 4 Hierarchical MAP estimation with unknown regularisation parameters Conclusion 5 M. Pereyra (MI — HWU) Bayesian mathematical imaging 17 / 53
Hypothesis testing Bayesian hypothesis test for specific image structures (e.g., lesions) H 0 ∶ The structure of interest is ABSENT in the true image H 1 ∶ The structure of interest is PRESENT in the true image The null hypothesis H 0 is rejected with significance α if P ( H 0 ∣ y ) ≤ α. Theorem (Repetti et al., 2018) Let S denote the region of R d associated with H 0 , containing all images without the structure of interest. Then S ∩ ̃ C α = ∅ � ⇒ P ( H 0 ∣ y ) ≤ α . If in addition S is convex, then checking S ∩ ̃ C α = ∅ is a convex problem x ∈ ̃ x , x ∈ R d ∥ ¯ x − x ∥ 2 x ∈ S . min s.t. ¯ C α , 2 ¯ M. Pereyra (MI — HWU) Bayesian mathematical imaging 18 / 53
Uncertainty quantification in MRI imaging x ∈ ̃ x ∈ S x MAP ˆ ¯ C 0 . 01 x ∈ ̃ ˆ x MAP (zoom) ¯ C 0 . 01 (zoom) x ∈ S (zoom) x = x, hence we fail to reject H 0 and conclude that MRI experiment: test images ¯ there is little evidence to support the observed structure. M. Pereyra (MI — HWU) Bayesian mathematical imaging 19 / 53
Uncertainty quantification in MRI imaging ˆ ¯ x ∈ C 0 . 01 x ∈ S 0 x MAP ˆ x MAP (zoom) ¯ x ∈ C 0 . 01 (zoom) x ∈ S 0 (zoom) x ≠ x, hence we reject H 0 and conclude that there is MRI experiment: test images ¯ significant evidence in favour of the observed structure. M. Pereyra (MI — HWU) Bayesian mathematical imaging 20 / 53
Recommend
More recommend
