scispace - formally typeset

Journal ArticleDOI

Convergence in Variance of Chebyshev Accelerated Gibbs Samplers

04 Feb 2014-SIAM Journal on Scientific Computing (Society for Industrial and Applied Mathematics)-Vol. 36, Iss: 1

TL;DR: An algorithm for the stochastic version of the second-order Chebyshev accelerated SSOR (symmetric successive overrelaxation) iteration is given and numerical examples of sampling from multivariate Gaussian distributions are provided to confirm that the desired convergence properties are achieved in finite precision.
Abstract: A stochastic version of a stationary linear iterative solver may be designed to converge in distribution to a probability distribution with a specified mean μ and covariance matrix A−1. A common example is Gibbs sampling applied to a multivariate Gaussian distribution which is a stochastic version of the Gauss-Seidel linear solver. The iteration operator that acts on the error in mean and covariance in the stochastic iteration is the same iteration operator that acts on the solution error in the linear solver, and thus both the stationary sampler and the stationary solver have the same error polynomial and geometric convergence rate. The polynomial acceleration techniques that are well known in numerical analysis for accelerating the linear solver may also be used to accelerate the stochastic iteration. We derive first-order and second-order Chebyshev polynomial acceleration for the stochastic iteration to accelerate convergence in the mean and covariance by mimicking the derivation for the linear solver. In particular, we show that the error polynomials are identical and hence so are the convergence rates. Thus, optimality of the Chebyshev accelerated solver implies optimality of the Chebyshev accelerated sampler. We give an algorithm for the stochastic version of the second-order Chebyshev accelerated SSOR (symmetric successive overrelaxation) iteration and provide numerical examples of sampling from multivariate Gaussian distributions to confirm that the desired convergence properties are achieved in finite precision.
Topics: Chebyshev iteration (70%), Chebyshev nodes (62%), Solver (61%), Chebyshev equation (60%), Chebyshev polynomials (59%)

Content maybe subject to copyright    Report

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
SIAM J. SCI. COMPUT.
c
2014 Society for Industrial and Applied Mathematics
Vol. 36, No. 1, pp. A124–A147
CONVERGENCE IN VARIANCE OF CHEBYSHEV ACCELERATED
GIBBS SAMPLERS
COLIN FOX
AND ALBERT PARKER
Abstract. A stoc hastic version of a stationary linear iterativ e solver may be designed to converge
in distribution to a probability distribution with a specified mean μ and covariance matrix A
1
.
A common example is Gibbs sampling applied to a multivariate Gaussian distribution which is a
stochastic version of the Gauss–Seidel linear solver. The iteration operator that acts on the error in
mean and covariance in the stochastic iteration is the same iteration operator that acts on the solution
error in the linear solver, and thus both the stationary sampler and the stationary solv er have the
same error polynomial and geometric convergence rate. The polynomial acceleration techniques that
are well known in numerical analysis for accelerating the linear solver may also be used to accelerate
the stochastic iteration. We derive first-order and second-order Chebyshev polynomial acceleration
for the stochastic iteration to accelerate convergence in the mean and covariance by mimicking the
derivation for the linear solver. In particular, we show that the error polynomials are identical and
hence so are the con vergence rates. Thus, optimality of the Chebyshev accelerated solver implies
optimality of the Chebyshev accelerated sampler. We give an algorithm for the stochastic version of
the second-order Cheb yshev accelerated SSOR (symmetric successive overrelaxation) iteration and
provide numerical examples of sampling from multivariate Gaussian distributions to confirm that
the desired convergence properties are achieved in finite precision.
Key words. Chebyshev polynomial acceleration, Gauss–Seidel, Gibbs sampling, geometric
convergence, linear solver, stochastic iteration, SSOR
AMS subject classifications. 65F10, 65B99, 62E17, 60G15, 60G60
DOI. 10.1137/120900940
1. Introduction. Iterations of the form
(1.1) x
l+1
= Gx
l
+ g, l =1, 2,...,
where G is a fixed iteration operator and g is a fixed vector, are commonplace in
numerical computation. For example, they occur in the stationary linear iterative
methods used to solve systems of linear equations [1, 10, 17, 23]. We often refer to the
associated algorithm as a solver. We consider these iterations, and also the related
stochastic iteration
(1.2) y
l+1
= Gy
l
+ g
l
,l=1, 2,...,
where now g
l
is a “noise” vector given by an independent draw from some fixed
probability distribution with finite variance. Just as the deterministic iteration (1.1)
can be designed to converge to the solution of a linear system that is too large or
complex to solve directly, the stochastic iteration (1.2) may be designed to converge
in distribution to a target distribution that is too high dimensional, or complex, to
sample from directly. Since the stochastic iteration may be used to generate samples
Submitted to the journal’s Methods and Algorithms for Scientific Computing section Decem-
ber 3, 2012; accepted for publication (in revised form) November 14, 2013; published electronically
February 4, 2014. This work was supported by the New Zealand Institute for Mathematics and its
Applications thematic programme on PDEs and Marsden contract UOO1015.
http://www.siam.org/journals/sisc/36-1/90094.html
Department of Physics, University of Otago, Dunedin, New Zealand (fox@physics.otago.ac.nz).
Cent er for Biofilm Engineering and Department of Mathematical Sciences, Montana State Uni-
versit y, Bozeman, MT 59715 (parker@math.montana.edu).
A124
Downloaded 03/07/14 to 153.90.118.169. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
CHEBYSHEV SAMPLING A125
from a desired target distribution, we often refer to the associated algorithm as a
sampler. An example is the conventional Gibbs sampling algorithm [21] applied to
sampling from a high-dimensional Gaussian distribution. In that case the iteration
operator G is identical to the iteration operator in the Gauss–Seidel iterative method
[5, 7].
Novel Gibbs samplers may be designed by considering matrix splittings other than
the Gauss–Seidel splitting [5]. Matrix splittings are considered further in section 2.
Interestingly, the deterministic and stochastic iterations converge under exactly the
same conditions, with a necessary and sufficient condition being that the spectral
radius of G be strictly less than 1, that is, ρ (G) < 1 [4, 26]. Convergence in both cases
is geometric, with the asymptotic average reduction factor given by ρ (G) (though this
is called the “convergence rate” in the statistics literature [19]).
A standard method of reducing the asymptotic average reduction factor is by
polynomial acceleration, particularly using Chebyshev polynomials [1, 6, 10, 23]. The
original formulation used a modified first-order iteration, as above, though the result-
ing algorithm is impractical due to numerical difficulties [1]. Practical implementa-
tions use a nonstationary second-order iteration that can give optimal reduction of
error at each iteration.
In this paper, we develop polynomial acceleration for the stochastic iteration.
In particular, we develop nonstationary first- and second-order iterations that give
optimal convergence in mean and variance to a desired target distribution. Since
convergence in mean is achieved by using exactly the linear iteration for solving a linear
system, polynomial acceleration of the mean is exactly as in the existing treatments.
Hence we focus on optimal convergence in variance that requires modification to the
noise term. Correspondingly, we focus throughout the development on sampling from
a target distribution that has zero mean and some finite covariance matrix, and hence
the noise distribution always has zero mean. Extension to target distributions with
nonzero mean is achieved simply by adding the deterministic iteration or, equivalently,
adding a fixed vector to the noise term.
We develop the sampling algorithms and demonstrate the equivalence to linear
solvers by investigating a sequence of linear iterative solvers, essentially following the
historical development in sophistication and speed, and show that exactly the same
ideas used to establish properties of the solver can be used to establish the equivalent
properties for a sampler. In particular, convergence of the solver implies convergence
of the sampler, and the convergence factors are identical, because they are given by
the same expression.
We follow the development and derivations of convergence, given in Axelsson [1],
for stationary and nonstationary (Chebyshev) first-order and second-order methods,
set out in sect. 5.2 (Stationary Iterative Methods) and sect. 5.3 (The Chebyshev It-
erative Method). We could have equally followed the excellent presentations of the
same methods in Golub and Van Loan [10] or Saad [23]. Our own work and compu-
tational implementation actually take a route that switches between the formalism
used in these three texts. By following here the route of a single exposition, we hope
to show how establishing convergence of the stochastic versions can be made very
straightforward.
The most straightforward application of the methods we develop is to sample from
a high-dimensional Gaussian distribution, defined by the mean vector μ and covariance
matrix A
1
. We present an example which shows the convergence of the Chebyshev
sampler in finite precision applied to a Gaussian Markov random field (GMRF) with
a known sparse precision matrix corresponding to a Mat´ern-class covariance function
Downloaded 03/07/14 to 153.90.118.169. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
A126 COLIN FOX AND ALBERT PARKER
[12, 15]. This example allows efficient numerical calculation since operation by A has
reduced numerical cost.
Although we focus on the Gaussian in our numerical example, the accelerated
algorithms we give are more generally applicable to any distribution where the focus
is on the mean as a “best” estimate and the covariance as a measure of uncertainties,
with higher moments not of primary concern. This is typical in inferential methods
applied to solving inverse problems or in the growing field of uncertainty quantifica-
tion, where the mean and variance of the distribution over parameters or predicted
quantities are the primary summary statistics of interest.
1.1. Some links between sampling from distributions and solving sys-
tems of equations. Consider a probability distribution with probability density
function π(x) and the two tasks of drawing x π (x distributed as π)andcom-
puting x =argmaxπ (or solving −∇ log π = 0). We use the notation x
i
=
(x
1
,x
2
,...,x
i1
,x
i+1
,...,x
n
)todenotealln 1componentsofx other than x
i
,
and π (x
i
|x
i
) to denote the univariate conditional distribution over x
i
conditioned
on the (fixed) value of all other components.
The classical Gibbs sampler or “stochastic relaxation” (also known as Glauber
dynamics and the local heat bath algorithm) for generating a sample from π is an
iterative algorithm in which one sweep consists of updating each component in se-
quence by drawing from the conditional distribution for the component with all other
components fixed at the most recent value, as in Algorithm 1. Repeating this sweep
indefinitely produces distributions over iterates that are guaranteed to converge (ge-
ometrically) to π under mild conditions [11, ref. 84], [19], though distributions with
nonconnected support for which Algorithm 1 fails are easy to find [19].
Algorithm 1: One sweep of the componentwise Gibbs sampler targeting π(x)
for i =1,...,n do
sample z π (x
i
|x
i
);
x
i
= z;
end
It is not hard to see a connection between the Gibbs sampler in Algorithm 1 and
the traditional Gauss–Seidel algorithm for maximizing π which consists of repeatedly
applying the sweep over componentwise solvers with all other components fixed at the
most recent value, as in Algorithm 2: Whereas the Gibbs sampler performs a compo-
nentwise conditional sampling, Gauss–Seidel performs componentwise optimization.
Algorithm 2: One sweep of Gauss–Seidel relaxation for maximizing π(x)
for i =1,...,n do
set z =argmax
x
i
π (x
i
|x
i
);
x
i
= z;
end
In statistical physics, distributions often arise with the form
π(x)=k(β)e
βH(x)
,
Downloaded 03/07/14 to 153.90.118.169. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
CHEBYSHEV SAMPLING A127
where H(x) is an energy function (the Hamiltonian), β is inversely proportional to
temperature, and k is a normalizing constant. It is often noted that a sampling
algorithm may be used to minimize H(x) in the zero temperature limit, i.e., by taking
the limit β →∞. Then sampling degenerates to optimization since the distribution
is localized at the mode. In particular, Algorithm 1 reduces to Algorithm 2.
In this paper, we exploit an equivalence that operates at finite β to show how the
minimizer (or solver) may be adapted to become a sampling algorithm. For example,
in the simplest case that β =1andH is quadratic, i.e.,
H(x)=
1
2
x
T
Ax b
T
x
for some symmetric positive definite (precision matrix) A, π is Gaussian and the
Gauss–Seidel minimizer of H becomes the Gibbs sampler for π when coordinatewise
minimization is replaced by coordinatewise conditional sampling. One sweep of the
Gibbs sampler may be written in the matrix form (1.2) with
G = M
1
N and g
l
= M
1
c
l
, where c
l
iid
N(0,D).
Here M = L+ D and N = L
T
is a splitting of the (symmetric) precision matrix A in
which L is the strictly lower triangular part of A and D is the diagonal of A [11]. This
is the same splitting used to write the Gauss–Seidel algorithm for solving Ax = b in
matrix form (1.1), with g = M
1
b. What makes this correspondence important is that
the convergence properties of the solver are inherited by the sampler (and vice versa),
which means that acceleration techniques developed for the solver may be applied
to the sampler. The main purpose of this paper is to establish the equivalence of
convergence in mean and covariance in the case of Chebyshev polynomial acceleration,
without the assumption of the target distribution being Gaussian.
2. Matrix splitting and iteration operators. Consider the splitting
(2.1) A = M N,
where A is a symmetric positive definite (SPD) matrix and M is invertible. For
example, for the Gauss–Seidel iteration, M is set to the lower triangular part of
A (including the diagonal). We will often consider the case where the splitting is
symmetric, which means that M is symmetric, and hence so is N. We will utilize the
family of iteration operators
(2.2) G
τ
=
I τM
1
A
parameterized by the relaxation parameter τ = 0. The natural iteration operator
induced by the splitting (2.1) is the case τ = 1, which we denote by G. The nonsta-
tionary iterative methods that we consider use a sequence of iteration operators with
parameters τ
l
, l =0, 1, 2,..., where l denotes iteration number. We will abbreviate
G
τ
l
by G
l
, where possible, to avoid subscripts on subscripts.
The iteration operator G
τ
may also be thought of as being induced by the splitting
A = M
τ
N
τ
,(2.3)
where
M
τ
=
1
τ
M and N
τ
= N +
1 τ
τ
M.
Downloaded 03/07/14 to 153.90.118.169. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
A128 COLIN FOX AND ALBERT PARKER
In the remainder of this section we list some lemmas about iteration operators
that we will use. Throughout the rest of the paper, proofs to lemmas and some
theorems have been deferred to the appendix.
Lemma 2.1. The itera tion operators G
τ
and G
κ
commute, that is, G
τ
G
κ
= G
κ
G
τ
for all τ, κ.
Lemma 2.2. For a symmetric splitting, G
τ
A
1
is symmetric.
The following lemma determines the variance of noise terms in sampling algo-
rithms.
Lemma 2.3. A
1
G
τ
A
1
G
T
τ
= M
1
τ
M
T
τ
+ N
τ
M
T
τ
.
3. First-order iterative methods. We first consider iterative solvers of the
equation
Ax
= b,
where A is a given SPD matrix, b is a given vector, and the solution we seek is denoted
by x
.
3.1. First-order stationary iterative solver. The first-order stationary iter-
ative solver uses the iteration
x
l+1
= x
l
τM
1
r
l
= G
τ
x
l
+ g
τ
,(3.1)
where r
l
= Ax
l
b for l =1, 2,..., the iteration operator G
τ
is given by (2.2)
and g
τ
= τM
1
b. In the remainder of this section, we derive the fixed point, error
polynomial, and average reduction factor for this iteration.
Lemma 3.1. The iteration in (3.1) has x
as its unique fixed point, i.e.,
(3.2) x
= G
τ
x
+ g
τ
Ax
= b.
Define the error at the lth iteration by
(3.3) e
l
= x
l
x
.
Subtract (3.2) from (3.1) to get the iteration for error
e
l+1
= x
l+1
x
= G
τ
x
l
+ g
τ
G
τ
x
g
τ
= G
τ
x
l
x
= G
τ
e
l
.
By recursion we prove the following theorem.
Theorem 3.2.
e
m
= G
m
τ
e
0
=
I τM
1
A
m
e
0
= P
m
M
1
A
e
0
,
where P
m
is the (simple) mth-order polynomial P
m
(λ)=(1 τλ)
m
.
Note that P
m
(0) = 1 and P
m
(1) = 0. The convergence and convergence rate
of the stationary iterative solver follow from Theorem 3.2.
Axelsson [1, p. 176] gives the optimal relaxation parameter
τ
opt
=
2
λ
1
+ λ
n
,
where λ
1
n
are the extreme (positive) eigenvalues of M
1
A, giving the average
reduction factor
(3.4) ρ
0
=
1 λ
1
n
1+λ
1
n
.
Note that this implies that the iterative solver (3.1) converges for some value of τ.To
be more precise, as long as M
1
A has all positive eigenvalues, then λ
1
n
(0, 1),
which means that ρ
0
(0, 1) and the iteration is guaranteed to converge.
Downloaded 03/07/14 to 153.90.118.169. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Citations
More filters

Journal ArticleDOI
Alessio Spantini1, Antti Solonen2, Tiangang Cui1, James Martin3  +2 moreInstitutions (4)
Abstract: In the Bayesian approach to inverse problems, data are often informative, relative to the prior, only on a low-dimensional subspace of the parameter space. Significant computational savings can be achieved by using this subspace to characterize and approximate the posterior distribution of the parameters. We first investigate approximation of the posterior covariance matrix as a low-rank update of the prior covariance matrix. We prove optimality of a particular update, based on the leading eigendirections of the matrix pencil defined by the Hessian of the negative log-likelihood and the prior precision, for a broad class of loss functions. This class includes the Forstner metric for symmetric positive definite matrices, as well as the Kullback--Leibler divergence and the Hellinger distance between the associated distributions. We also propose two fast approximations of the posterior mean and prove their optimality with respect to a weighted Bayes risk under squared-error loss. These approximations are dep...

104 citations


Proceedings Article
05 Dec 2013
TL;DR: It is shown that if the Gaussian precision matrix is generalized diagonally dominant, then any Hogwild Gibbs sampler, with any update schedule or allocation of variables to processors, yields a stable sampling process with the correct sample mean.
Abstract: Sampling inference methods are computationally difficult to scale for many models in part because global dependencies can reduce opportunities for parallel computation. Without strict conditional independence structure among variables, standard Gibbs sampling theory requires sample updates to be performed sequentially, even if dependence between most variables is not strong. Empirical work has shown that some models can be sampled effectively by going "Hogwild" and simply running Gibbs updates in parallel with only periodic global communication, but the successes and limitations of such a strategy are not well understood. As a step towards such an understanding, we study the Hogwild Gibbs sampling strategy in the context of Gaussian distributions. We develop a framework which provides convergence conditions and error bounds along with simple proofs and connections to methods in numerical linear algebra. In particular, we show that if the Gaussian precision matrix is generalized diagonally dominant, then any Hogwild Gibbs sampler, with any update schedule or allocation of variables to processors, yields a stable sampling process with the correct sample mean.

62 citations



Journal ArticleDOI
TL;DR: Standard Gibbs sampling applied to a multivariate normal distribution with a specified precision matrix is equivalent in fundamental ways to the Gauss-Seidel iterative solution of linear equations in the precision matrix, giving easy access to mature results in numerical linear algebra.
Abstract: Standard Gibbs sampling applied to a multivariate normal distribution with a specified precision matrix is equivalent in fundamental ways to the Gauss–Seidel iterative solution of linear equations in the precision matrix. Specifically, the iteration operators, the conditions under which convergence occurs, and geometric convergence factors (and rates) are identical. These results hold for arbitrary matrix splittings from classical iterative methods in numerical linear algebra giving easy access to mature results in that field, including existing convergence results for antithetic-variable Gibbs sampling, REGS sampling, and generalizations. Hence, efficient deterministic stationary relaxation schemes lead to efficient generalizations of Gibbs sampling. The technique of polynomial acceleration that significantly improves the convergence rate of an iterative solver derived from a symmetric matrix splitting may be applied to accelerate the equivalent generalized Gibbs sampler. Identicality of error polynomials guarantees convergence of the inhomogeneous Markov chain, while equality of convergence factors ensures that the optimal solver leads to the optimal sampler. Numerical examples are presented, including a Chebyshev accelerated SSOR Gibbs sampler applied to a stylized demonstration of low-level Bayesian image reconstruction in a large 3-dimensional linear inverse problem.

19 citations


Cites background or methods from "Convergence in Variance of Chebyshe..."

  • ...The Chebyshev accelerated SSOR solver and corresponding Chebyshev accelerated SSOR sampler (Fox and Parker, 2014) are depicted in panels C and D of Figure 1....

    [...]

  • ...But even sooner, after k∗∗ = k∗/2 iterations, the Chebyshev error reduction for the variance is predicted to be smaller than ε (Fox and Parker [19])....

    [...]

  • ...Using Theorem 5, we derived the Chebyshev accelerated SSOR sampler (Fox and Parker, 2014) by iteratively updating parameters via (13) and then generating a sampler via (17)....

    [...]

  • ...But even sooner, after k∗∗ = k∗/2 iterations, the Chebyshev error reduction for the variance is predicted to be smaller than ε (Fox and Parker, 2014). imsart-bj ver....

    [...]

  • ...Fox and Parker (2014) considered point-wise convergence of the mean and variance of a Gibbs SSOR sampler accelerated by Chebyshev polynomials....

    [...]


Dissertation
01 Jan 2017
TL;DR: A rigorous mathematical characterization of lowdimensional structures that enable efficient Bayesian inference in high-dimensional and continuous parameter spaces is provided and this characterization is exploited to devise new structure-exploiting and computationally efficient inference algorithms.
Abstract: The Bayesian approach to inference characterizes model parameters and predictions through the exploration of their posterior distributions, i.e., their distributions conditioned on available data. The Bayesian paradigm provides a flexible, principled framework for quantifying uncertainty, wherein heterogeneous and incomplete sources of information (e.g., prior knowledge, noisy observations, imperfect models) can be properly rationalized. Yet a major obstacle to deploying Bayesian inference in realistic applications is computational: characterizing the associated high-dimensional and non-Gaussian posterior distributions remains a challenging task. While the Bayesian formulation is quite general, essential features of a statistical model can bring additional structure to the Bayesian update. For instance, the prior distribution often encodes some kind of regularity in the parameters; observations might be sparse and corrupted by noise; observations might also be indirect, related to the parameters by a forward operator that filters out some information; the posterior distribution might satisfy conditional independence assumptions that reflect local probabilistic interactions; and in some cases we might be uninterested in the posterior distribution per se, but rather in specific prediction goals. In this thesis we: (1) provide a rigorous mathematical characterization of lowdimensional structures that enable efficient Bayesian inference in high-dimensional and continuous parameter spaces; and (2) exploit this characterization to devise new structure-exploiting and computationally efficient inference algorithms. Our contributions encompass multiple related topics. First we characterize optimal low-rank approximations of linear–Gaussian Bayesian inverse problems, and of their goal-oriented extensions. Then we turn to inference in the nonlinear non-Gaussian setting—analyzing the sparsity, decomposability, and low-rank structure of deterministic couplings between distributions. These couplings facilitate efficient computation of posterior expectations in generically non-Gaussian settings. Based on this analysis, we introduce a number of approaches for representing non-Gaussian Markov random fields and for exploiting their conditional independence structure in computation by means of sparse nonlinear transport maps. We also develop new variational

14 citations


Cites background from "Convergence in Variance of Chebyshe..."

  • ...Alternative techniques for sampling from high-dimensional Gaussian distributions can be found, for instance, in [211, 100]....

    [...]


References
More filters

Book
01 Jan 1983

34,706 citations


"Convergence in Variance of Chebyshe..." refers methods in this paper

  • ...A standard method of reducing the asymptotic average reduction factor is by polynomial acceleration, particularly using Chebyshev polynomials [1, 6, 10, 23]....

    [...]

  • ...Analogous to the SSOR solver algorithms in [10, 23], the Chebyshev sampler implements sequential forward and backward sweeps of an SOR sampler [7, 20] (i....

    [...]

  • ...Since a symmetric splitting is required, Chebyshev acceleration in a linear solver is commonly implemented with a symmetric successive overrelaxation (SSOR) splitting A = MSSOR−NSSOR, with algorithms to be found, for example, in [10, 23]....

    [...]

  • ...For example, they occur in the stationary linear iterative methods used to solve systems of linear equations [1, 10, 17, 23]....

    [...]

  • ...We could have equally followed the excellent presentations of the same methods in Golub and Van Loan [10] or Saad [23]....

    [...]


Journal ArticleDOI
Stuart Geman1, Donald Geman2Institutions (2)
TL;DR: The analogy between images and statistical mechanics systems is made and the analogous operation under the posterior distribution yields the maximum a posteriori (MAP) estimate of the image given the degraded observations, creating a highly parallel ``relaxation'' algorithm for MAP estimation.
Abstract: We make an analogy between images and statistical mechanics systems. Pixel gray levels and the presence and orientation of edges are viewed as states of atoms or molecules in a lattice-like physical system. The assignment of an energy function in the physical system determines its Gibbs distribution. Because of the Gibbs distribution, Markov random field (MRF) equivalence, this assignment also determines an MRF image model. The energy function is a more convenient and natural mechanism for embodying picture attributes than are the local characteristics of the MRF. For a range of degradation mechanisms, including blurring, nonlinear deformations, and multiplicative or additive noise, the posterior distribution is an MRF with a structure akin to the image model. By the analogy, the posterior distribution defines another (imaginary) physical system. Gradual temperature reduction in the physical system isolates low energy states (``annealing''), or what is the same thing, the most probable states under the Gibbs distribution. The analogous operation under the posterior distribution yields the maximum a posteriori (MAP) estimate of the image given the degraded observations. The result is a highly parallel ``relaxation'' algorithm for MAP estimation. We establish convergence properties of the algorithm and we experiment with some simple pictures, for which good restorations are obtained at low signal-to-noise ratios.

18,328 citations


Book
01 Apr 2003
TL;DR: This chapter discusses methods related to the normal equations of linear algebra, and some of the techniques used in this chapter were derived from previous chapters of this book.
Abstract: Preface 1. Background in linear algebra 2. Discretization of partial differential equations 3. Sparse matrices 4. Basic iterative methods 5. Projection methods 6. Krylov subspace methods Part I 7. Krylov subspace methods Part II 8. Methods related to the normal equations 9. Preconditioned iterations 10. Preconditioning techniques 11. Parallel implementations 12. Parallel preconditioners 13. Multigrid methods 14. Domain decomposition methods Bibliography Index.

12,575 citations


"Convergence in Variance of Chebyshe..." refers methods in this paper

  • ...A standard method of reducing the asymptotic average reduction factor is by polynomial acceleration, particularly using Chebyshev polynomials [1, 6, 10, 23]....

    [...]

  • ...Analogous to the SSOR solver algorithms in [10, 23], the Chebyshev sampler implements sequential forward and backward sweeps of an SOR sampler [7, 20] (i....

    [...]

  • ...Since a symmetric splitting is required, Chebyshev acceleration in a linear solver is commonly implemented with a symmetric successive overrelaxation (SSOR) splitting A = MSSOR−NSSOR, with algorithms to be found, for example, in [10, 23]....

    [...]

  • ...For example, they occur in the stationary linear iterative methods used to solve systems of linear equations [1, 10, 17, 23]....

    [...]

  • ...We could have equally followed the excellent presentations of the same methods in Golub and Van Loan [10] or Saad [23]....

    [...]


Book
Gene H. Golub1, Charles Van Loan2Institutions (2)
01 Nov 1996

8,491 citations


Book
Christian P. Robert1, George Casella1Institutions (1)
01 Jan 1999
TL;DR: This new edition contains five completely new chapters covering new developments and has sold 4300 copies worldwide of the first edition (1999).
Abstract: We have sold 4300 copies worldwide of the first edition (1999). This new edition contains five completely new chapters covering new developments.

6,854 citations


"Convergence in Variance of Chebyshe..." refers background in this paper

  • ...84], [19], though distributions with nonconnected support for which Algorithm 1 fails are easy to find [19]....

    [...]

  • ...Convergence in both cases is geometric, with the asymptotic average reduction factor given by ρ (G) (though this is called the “convergence rate” in the statistics literature [19])....

    [...]