krylov_sampling.py 2.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2018 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.

from numpy import sqrt
20
from numpy.random import randn
21
22


Martin Reinecke's avatar
Martin Reinecke committed
23
24
def generate_krylov_samples(D_inv, S, j=None,  N_samps=1, N_iter=10,
                            name=None):
25
    """
Martin Reinecke's avatar
Martin Reinecke committed
26
    Generates inverse samples from a curvature D.
27
28
29
30
31
32
33
34
    This algorithm iteratively generates samples from
    a curvature D by applying conjugate gradient steps
    and resampling the curvature in search direction.

    Parameters
    ----------
    D_inv : EndomorphicOperator
        The curvature which will be the inverse of the covarianc
35
        of the generated samples 
36
37
38
39
40
41
42
    S : EndomorphicOperator (from which one can sample)
        A prior covariance operator which is used to generate prior
        samples that are then iteratively updated
    j : Field, optional
        A Field to which the inverse of D_inv is applied. The solution
        of this matrix inversion problem is a side product of generating
        the samples.
43
        If not supplied, it is sampled from the inverse prior.
44
45
46
47
48
49
50
51
52
    N_samps : Int, optional
        How many samples to generate. Default: 1
    N_iter : Int, optional
        How many iterations of the conjugate gradient to run. Default: 10

    Returns
    -------
    (solution, samples) : A tuple of a field 'solution' and a list of fields
        'samples'. The first entry of the tuple is the solution x to
53
            D_inv(x) = j
54
55
        and the second entry are a list of samples from D_inv.inverse
    """
Martin Reinecke's avatar
Martin Reinecke committed
56
    j = S.draw_sample(from_inverse=True) if j is None else j
57
    x = j*0
58
59
60
    r = j.copy()
    p = r.copy()
    d = p.vdot(D_inv(p))
Martin Reinecke's avatar
Martin Reinecke committed
61
62
63
    y = [S.draw_sample() for _ in range(N_samps)]
    for k in range(1, 1+N_iter):
        gamma = r.vdot(r)/d
64
65
        if gamma == 0.:
            break
Martin Reinecke's avatar
Martin Reinecke committed
66
        x += gamma*p
67
68
69
70
        for i in range(N_samps):
            y[i] -= p.vdot(D_inv(y[i])) * p / d
            y[i] += randn() / sqrt(d) * p
        r_new = r - gamma * D_inv(p)
Martin Reinecke's avatar
Martin Reinecke committed
71
        beta = r_new.vdot(r_new) / r.vdot(r)
72
73
74
75
76
        r = r_new
        p = r + beta * p
        d = p.vdot(D_inv(p))
        if d == 0.:
            break
77
78
        if name is not None:
            print('{}: Iteration #{}'.format(name, k))
79
    return x, y