%matplotlib inline
%config InlineBackend.figure_format = 'svg'

import numpy as np
import scipy
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.gridspec as gridspec
import seaborn as sns

# Set matplotlib and seaborn plotting style
sns.set_style('darkgrid')
np.random.seed(42)

1. A kernel example

def exp_quadratic(xa, xb):
    """Exponentiated quadratic  with σ=1"""
    # L2 distance (Squared Euclidian)
    sq_norm = -0.5 * scipy.spatial.distance.cdist(xa, xb, 'sqeuclidean')
    return np.exp(sq_norm)  
X = np.expand_dims(np.linspace(*xlim, 25), 1)
Σ = exp_quadratic(X, X)
plt.imshow(Σ, cmap=cm.YlGnBu);
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"> 2022-03-28T15:57:12.744982 image/svg+xml Matplotlib v3.5.1, https://matplotlib.org/
zero = np.array([[0]])
Σ0 = exp_quadratic(X, zero)
plt.plot(X[:,0], Σ0[:,0]);
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"> 2022-03-28T15:57:14.116661 image/svg+xml Matplotlib v3.5.1, https://matplotlib.org/

2. Sampling from prior

n_samples = 100
n_funcs = 8

X = np.expand_dims(np.linspace(-4,4, n_samples), 1)
Σ = exp_quadratic(X, X)

ys = np.random.multivariate_normal(mean=np.zeros(n_samples), cov=Σ, size=n_funcs)
for i in range(n_funcs):
    plt.plot(X, ys[i], linestyle='-', marker='o', markersize=3)
plt.xlabel('$x$', fontsize=13)
plt.ylabel('$y = f(x)$', fontsize=13)
plt.title((
    f'{n_funcs} different function realizations at {n_samples} points\n'
    'sampled from a Gaussian process with exponentiated quadratic kernel'))
plt.xlim([-4, 4])
plt.show()
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"> 2022-03-28T16:01:41.515158 image/svg+xml Matplotlib v3.5.1, https://matplotlib.org/
exponentiated_quadratic = exp_quadratic
A = np.array([[1,-2j],[2j,5]])
A, A.shape
(array([[ 1.+0.j, -0.-2.j],
        [ 0.+2.j,  5.+0.j]]),
 (2, 2))
Cholesky decomposition in numpy
L = np.linalg.cholesky(A)
np.dot(L, L.T.conj())
array([[1.+0.j, 0.-2.j],
       [0.+2.j, 5.+0.j]])
A = [[1,-2j],[2j,5]] # what happens if A is only array_like?
np.linalg.cholesky(A) # an ndarray object is returned
array([[1.+0.j, 0.+0.j],
       [0.+2.j, 1.+0.j]])
np.linalg.cholesky(np.matrix(A))
matrix([[1.+0.j, 0.+0.j],
        [0.+2.j, 1.+0.j]])
def GP(X1, y1, X2, kernel_func):
    cov11 = kernel_func(X1, X1)
    cov12 = kernel_func(X1, X2)
    solved = scipy.linalg.solve(cov11, cov12, assume_a='pos').T
    
    mu2 = solved @ y1
    cov22 = kernel_func(X2, X2)
    cov2 = cov22 - (solved @ cov12)
    return mu2, cov2
    
TODO: Convert K-1 by using cholesky

https://jaketae.github.io/study/gaussian-process/

def GP2(X1, y1, X2, kernel_func):
    K11 = kernel_func(X1, X1)
    K12 = kernel_func(X1, X2)
    K22 = kernel_func(X2, X2)
    
    #L = np.linalg.cholesky(K11)    
    mu2 = K12.T.dot(np.linalg.inv(K11)).dot(y1)
    cov2 = K22 - K12.T.dot(np.linalg.inv(K11)).dot(K12)
    return mu2, cov2
ny = 10 # Number of functions
domain = (-6, 6)
domain[0]+2, domain[1]-2, (n1, 1)
X1 = np.random.uniform(domain[0]+2, domain[1]-2, size=(n1, 1))
X1.shape
(50, 1)
%%prun
f_sin = lambda x: (np.sin(x)).flatten()

n1 = 40 # Train points
n2 = 75 # Test points
ny = 5 # Number of functions
domain = (-6, 6)

X1 = np.random.uniform(domain[0]+2, domain[1]-2, size=(n1, 1))
y1 = f_sin(X1)

X2 = np.linspace(domain[0], domain[1], n2).reshape(-1, 1)

# mu2, cov2 = GP(X1, y1, X2, exp_quadratic)
mu2, cov2 = GP2(X1, y1, X2, exp_quadratic)
sigma2 = np.sqrt(np.diag(cov2))

y2 = np.random.multivariate_normal(mean=mu2, cov=cov2, size=ny)
 
<string>:16: RuntimeWarning: invalid value encountered in sqrt
<string>:18: RuntimeWarning: covariance is not positive-semidefinite.
         296 function calls (287 primitive calls) in 0.044 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        2    0.040    0.020    0.040    0.020 linalg.py:476(inv)
        1    0.001    0.001    0.001    0.001 linalg.py:1482(svd)
        3    0.001    0.000    0.001    0.000 3828527866.py:1(exp_quadratic)
        1    0.000    0.000    0.002    0.002 {method 'multivariate_normal' of 'numpy.random.mtrand.RandomState' objects}
        3    0.000    0.000    0.000    0.000 {built-in method scipy.spatial._distance_pybind.cdist_sqeuclidean}
        1    0.000    0.000    0.044    0.044 {built-in method builtins.exec}
     17/8    0.000    0.000    0.042    0.005 {built-in method numpy.core._multiarray_umath.implement_array_function}
        1    0.000    0.000    0.000    0.000 {method 'flatten' of 'numpy.ndarray' objects}
        4    0.000    0.000    0.000    0.000 {method 'dot' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.044    0.044 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 function_base.py:23(linspace)
        3    0.000    0.000    0.000    0.000 socket.py:480(send)
        1    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(linspace)
        1    0.000    0.000    0.000    0.000 {method 'uniform' of 'numpy.random.mtrand.RandomState' objects}
        1    0.000    0.000    0.041    0.041 624826088.py:1(GP2)
        1    0.000    0.000    0.000    0.000 numeric.py:2344(within_tol)
        4    0.000    0.000    0.000    0.000 {method 'reduce' of 'numpy.ufunc' objects}
        2    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(result_type)
        3    0.000    0.000    0.000    0.000 distance.py:2616(cdist)
        6    0.000    0.000    0.000    0.000 {method 'astype' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.000    0.000 numeric.py:2264(isclose)
       16    0.000    0.000    0.000    0.000 {built-in method numpy.array}
        4    0.000    0.000    0.000    0.000 fromnumeric.py:70(_wrapreduction)
        2    0.000    0.000    0.000    0.000 warnings.py:35(_formatwarnmsg_impl)
        3    0.000    0.000    0.000    0.000 iostream.py:208(schedule)
        2    0.000    0.000    0.000    0.000 iostream.py:502(write)
        3    0.000    0.000    0.000    0.000 linalg.py:135(_commonType)
        1    0.000    0.000    0.000    0.000 <string>:1(<lambda>)
        2    0.000    0.000    0.000    0.000 _ufunc_config.py:32(seterr)
        2    0.000    0.000    0.000    0.000 warnings.py:403(__init__)
        2    0.000    0.000    0.000    0.000 {built-in method builtins.abs}
        3    0.000    0.000    0.000    0.000 linalg.py:107(_makearray)
        1    0.000    0.000    0.000    0.000 fromnumeric.py:1513(diagonal)
        2    0.000    0.000    0.000    0.000 warnings.py:20(_showwarnmsg_impl)
        2    0.000    0.000    0.000    0.000 linecache.py:82(updatecache)
        2    0.000    0.000    0.000    0.000 _ufunc_config.py:132(geterr)
        1    0.000    0.000    0.000    0.000 twodim_base.py:229(diag)
        1    0.000    0.000    0.000    0.000 {built-in method numpy.arange}
        3    0.000    0.000    0.000    0.000 fromnumeric.py:2355(all)
        2    0.000    0.000    0.000    0.000 iostream.py:420(_is_master_process)
        3    0.000    0.000    0.000    0.000 threading.py:1071(is_alive)
        9    0.000    0.000    0.000    0.000 _asarray.py:23(asarray)
        3    0.000    0.000    0.000    0.000 linalg.py:102(get_linalg_error_extobj)
        6    0.000    0.000    0.000    0.000 _asarray.py:110(asanyarray)
        2    0.000    0.000    0.040    0.020 <__array_function__ internals>:2(inv)
        1    0.000    0.000    0.000    0.000 <frozen importlib._bootstrap>:1017(_handle_fromlist)
        3    0.000    0.000    0.000    0.000 threading.py:1017(_wait_for_tstate_lock)
        1    0.000    0.000    0.000    0.000 {method 'any' of 'numpy.generic' objects}
       12    0.000    0.000    0.000    0.000 {built-in method builtins.len}
        2    0.000    0.000    0.000    0.000 numerictypes.py:285(issubclass_)
        4    0.000    0.000    0.000    0.000 fromnumeric.py:71(<dictcomp>)
        2    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(dot)
        1    0.000    0.000    0.000    0.000 numeric.py:2186(allclose)
        3    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(all)
        3    0.000    0.000    0.000    0.000 {method 'acquire' of '_thread.lock' objects}
        2    0.000    0.000    0.000    0.000 linecache.py:15(getline)
        4    0.000    0.000    0.000    0.000 linalg.py:125(_realType)
        3    0.000    0.000    0.000    0.000 linalg.py:193(_assert_stacked_2d)
       12    0.000    0.000    0.000    0.000 {built-in method builtins.issubclass}
        2    0.000    0.000    0.000    0.000 {method 'reshape' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.000    0.000 numerictypes.py:359(issubdtype)
        1    0.000    0.000    0.000    0.000 fromnumeric.py:2256(any)
        7    0.000    0.000    0.000    0.000 {method 'get' of 'dict' objects}
        1    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(allclose)
        3    0.000    0.000    0.000    0.000 iostream.py:97(_event_pipe)
        2    0.000    0.000    0.000    0.000 warnings.py:96(_showwarnmsg)
        1    0.000    0.000    0.001    0.001 <__array_function__ internals>:2(svd)
        8    0.000    0.000    0.000    0.000 {built-in method builtins.isinstance}
        2    0.000    0.000    0.000    0.000 linecache.py:37(getlines)
        1    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(diag)
        2    0.000    0.000    0.000    0.000 {built-in method posix.getpid}
        6    0.000    0.000    0.000    0.000 linalg.py:112(isComplexType)
        4    0.000    0.000    0.000    0.000 {built-in method builtins.getattr}
        1    0.000    0.000    0.000    0.000 {method 'diagonal' of 'numpy.ndarray' objects}
        2    0.000    0.000    0.000    0.000 linalg.py:199(_assert_stacked_square)
        2    0.000    0.000    0.000    0.000 {built-in method numpy.seterrobj}
        1    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(diagonal)
        1    0.000    0.000    0.000    0.000 _methods.py:53(_any)
        2    0.000    0.000    0.000    0.000 iostream.py:439(_schedule_flush)
        2    0.000    0.000    0.000    0.000 warnings.py:117(_formatwarnmsg)
        1    0.000    0.000    0.000    0.000 numeric.py:1865(isscalar)
        1    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(any)
        1    0.000    0.000    0.000    0.000 _ufunc_config.py:433(__enter__)
        2    0.000    0.000    0.000    0.000 {method 'startswith' of 'str' objects}
        1    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(ndim)
        1    0.000    0.000    0.000    0.000 _ufunc_config.py:438(__exit__)
        4    0.000    0.000    0.000    0.000 {method '__array_prepare__' of 'numpy.ndarray' objects}
        4    0.000    0.000    0.000    0.000 {built-in method numpy.geterrobj}
        1    0.000    0.000    0.000    0.000 <__array_function__ internals>:2(isclose)
        1    0.000    0.000    0.000    0.000 fromnumeric.py:3106(ndim)
        1    0.000    0.000    0.000    0.000 _ufunc_config.py:429(__init__)
        3    0.000    0.000    0.000    0.000 {method 'lower' of 'str' objects}
        4    0.000    0.000    0.000    0.000 {method 'items' of 'dict' objects}
        3    0.000    0.000    0.000    0.000 threading.py:513(is_set)
        1    0.000    0.000    0.000    0.000 {built-in method builtins.hasattr}
        2    0.000    0.000    0.000    0.000 {method 'endswith' of 'str' objects}
        2    0.000    0.000    0.000    0.000 linalg.py:472(_unary_dispatcher)
        3    0.000    0.000    0.000    0.000 {method 'append' of 'collections.deque' objects}
        3    0.000    0.000    0.000    0.000 fromnumeric.py:2350(_all_dispatcher)
        2    0.000    0.000    0.000    0.000 multiarray.py:716(dot)
        2    0.000    0.000    0.000    0.000 multiarray.py:644(result_type)
        1    0.000    0.000    0.000    0.000 numeric.py:2182(_allclose_dispatcher)
        3    0.000    0.000    0.000    0.000 {built-in method builtins.callable}
        1    0.000    0.000    0.000    0.000 function_base.py:18(_linspace_dispatcher)
        1    0.000    0.000    0.000    0.000 linalg.py:1478(_svd_dispatcher)
        1    0.000    0.000    0.000    0.000 {built-in method _operator.index}
        1    0.000    0.000    0.000    0.000 numeric.py:2260(_isclose_dispatcher)
        1    0.000    0.000    0.000    0.000 fromnumeric.py:3102(_ndim_dispatcher)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        1    0.000    0.000    0.000    0.000 fromnumeric.py:2251(_any_dispatcher)
        1    0.000    0.000    0.000    0.000 fromnumeric.py:1509(_diagonal_dispatcher)
        1    0.000    0.000    0.000    0.000 twodim_base.py:225(_diag_dispatcher)
fig, (ax1, ax2) = plt.subplots(
    nrows=2, ncols=1, figsize=(6, 6))

# Plot the distribution of the function (mean, covariance)
ax1.plot(X2, f_sin(X2), 'b--', label='$sin(x)$')
ax1.fill_between(X2.flat, mu2-2*sigma2, mu2+2*sigma2, color='red', 
                 alpha=0.15, label='$2 \sigma_{2|1}$')
ax1.plot(X2, mu2, 'r-', lw=2, label='$\mu_{2|1}$')
ax1.plot(X1, y1, 'ko', linewidth=2, label='$(x_1, y_1)$')

# Plot some samples from this function
ax2.plot(X2, y2.T, '-')
ax2.set_xlabel('$x$', fontsize=13)
ax2.set_ylabel('$y$', fontsize=13)
ax2.set_title('5 different function realizations from posterior')
ax1.axis([domain[0], domain[1], -3, 3])
ax2.set_xlim([-6, 6])
plt.tight_layout()
plt.show()
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"> 2022-03-28T17:15:07.852920 image/svg+xml Matplotlib v3.5.1, https://matplotlib.org/