on 07-Jan-2018 (Sun)

In linear algebra, the Cholesky decomposition or Cholesky factorization is a decomposition of a Hermitian, positive-definite matrix into the product of a lower triangular matrix and its conjugate transpose, which is useful e.g. for efficient numerical solutions and Monte Carlo simulations.

Hermitian matrices can be understood as the complex extension of real symmetric matrices.

Kalman filtering, also known as linear quadratic estimation (LQE), is an algorithm that uses a series of measurements observed over time, containing statistical noise and other inaccuracies, and produces estimates of unknown variables that tend to be more accurate than those based on a single measurement alone, by estimating a joint probability distribution over the variables for each timeframe.

"The differentiation of computational graphs that used the Cholesky decomposition required a new Op, which we contributed." (Matthews et al 2017:3)

"We also contributed code that enabled GPU solving of matrix triangular systems, which in some cases is the bottleneck for approximate inference." (Matthews et al 2017:3)

"Note that all the MCMC based inference methods support a Bayesian prior on the hyperparameters, whereas all other methods assume a point estimate. Our main MCMC method is Hamiltonian Monte Carlo (Neal, 2010)." (Matthews et al 2017:4)

"The object-oriented paradigm, whilst very natural for the Python layer of the code, has a dierent emphasis to that of a computational graph which is arguably closer to a functional concept. The largely functional computational graph, and a object-oriented interface need to live cleanly together in GP ow. This is achieved through the Param class that allows parameters to be both properties of an ob ject that can be manipulated as such and Variables in a TensorFlow graph that can be optimized." (Matthews et al 2017:4)

Since the objective function is unknown, the Bayesian strategy (of optimisation) is to treat it as a random function and place a prior over it. The prior captures our beliefs about the behaviour of the function. After gathering the function evaluations, which are treated as data, the prior is updated to form the posterior distribution over the objective function. The posterior distribution, in turn, is used to construct an acquisition function (often also referred to as infill sampling criteria) that determines what the next query point should be.

The Cholesky decomposition of a Hermitian positive-definite matrix A is a decomposition of the form

where L is a lower triangular matrix with real and positive diagonal entries, and L* denotes the conjugate transpose of L.

Every Hermitian positive-definite matrix (and thus also every real-valued symmetric positive-definite matrix) has a unique Cholesky decomposition.^[2] If the matrix A is Hermitian and positive semi-definite, then it still has a decomposition of the form A = LL* if the diagonal entries of L are allowed to be zero.^[3]

When A has real entries, L has real entries as well, and the factorization may be written A = LL^T.^[4]

The Cholesky decomposition is unique when A is positive definite; there is only one lower triangular matrix L with strictly positive diagonal entries such that A = LL*. However, the decomposition need not be unique when A is positive semidefinite.

The converse holds trivially: if A can be written as LL* for some invertible L, lower triangular or otherwise, then A is Hermitian and positive definite.

In numerical analysis and linear algebra, LU decomposition (where 'LU' stands for 'lower upper', and also called LU factorization) factors a matrix as the product of a lower triangular matrix and an upper triangular matrix.

Computers usually solve square systems of linear equations using the LU decomposition, and it is also a key step when inverting a matrix, or computing the determinant of a matrix.

Expecatation propagation aims to minimise \(KL(p(x,\theta\mid y,m)\parallel q(x,\theta\mid m))\). The problem with this minimisation, is that \( p(x,\theta\mid y,m)\) is unknown, so the first thing we do is to force some structure on this distribution, by making assumptions. For many models with IID observation, the posterior can be factorised as
$$ p(\theta \mid D) = \frac{1}{p(D)} p(\theta) \prod_{i=1}^{N} p(x^{(i)} \mid \theta) $$
which is a product of prior and likelihood, normalised by the data likelihood. Ignoring the constant, we can rearrange the posterior as
$$ p(\theta \mid D) \propto \prod_{i=0}^{N} f_i(\theta) $$
with \(f_0(\theta) = p(\theta) \) and \( f_i(\theta) = p(x^{(i)} \mid \theta) .\) And now our goal is to find a factorable
$$q(\theta) = \prod_{i=0}^{N} \tilde{f_i}(\theta) $$
that is close to our posterior.
EP manages this with a Gibbs sampling like procedure by
$$ \min_{\tilde{f_i}(\theta)} KL \left( f_i(\theta) \prod_{j \neq i}^{N} \tilde{f_j}(\theta) \parallel \tilde{f_i}(\theta) \prod_{j \neq i}^{N} \tilde{f_j}(\theta) \right) $$
and updating one factor at a time. The complete proceduce is as follows:

1. ​​​​Initialise $f_0(\theta) \cdots f_N(\theta) $; Initialise $\tilde{f_0}(\theta) = f_0(\theta) $ and $\tilde{f_i}(\theta) = 1 $ for $i > 0$
2. Compute $q(\theta) = \prod_{i=0}^{N} \tilde{f_i}(\theta) $
3. for $i = 0 \cdots N $ do
   1. Deletion: delete current variable from joint distribution: $q_{\i}(\theta) \leftarrow \frac{q(\theta)}{\tilde{f_i}(\theta)} = \prod_{j \neq i}^{N} \tilde{f_j}(\theta) $
   2. Projection: update current variable by minimising the KL divergence $\min_{\tilde{f_i}(\theta)} KL \left( f_i(\theta) q_{\\i}(\theta) \parallel \tilde{f_i}(\theta) q_{\\i}(\theta) \right) $
   3. Inclusion: update $q(\theta)$ using the current variable distribution $q(\theta ) \leftarrow \tilde{f_i^{new}} (\theta) q_{\\i}(\theta)
      

EP finds approximations to a probability distribution. It uses an iterative approach that leverages the factorization structure of the target distribution.

VB generally requiring more iterations, however sometimes each iteration is much cheaper than an EP iteration.

Variational Bayes can be seen as an extension of the EM (expectation-maximization) algorithm from maximum a posteriori estimation (MAP estimation) of the single most probable value of each parameter to fully Bayesian estimation which computes (an approximation to) the entire posterior distribution of the parameters and latent variables.