Inverse Covariance Kernels in 1D

cov_inv_1D.py

import numpy as np
import matplotlib
matplotlib.rc('text', usetex=True)
import matplotlib.pyplot as plt

# Properties of kernel
a = 0.     # Coeff. on Gaussian term (exponential term gets 1-a)
ell = 2.    # Exponential correlation length
sigma = 1.  # Gaussian std. dev. parameter

# 1D pixel array
j = np.arange(-10,11,1)

# Kronecker Delta
d = np.zeros(j.size)
d[-1] = 1.

# Kernel
kern_gauss = np.exp(-0.5 * (np.abs(j)/sigma)**2.)
kern_exp = np.exp(-np.abs(j/ell))
kern = a*kern_gauss + (1.-a)*kern_exp

# Inverse Kernel from Fourier-space formula: K^-1 = d/K
kern_inv = np.real(np.fft.ifft(np.fft.fft(d) / np.fft.fft(kern)))

# Plot kernel and its inverse
fig = plt.figure(figsize=(5,5), dpi=150)

ax = fig.add_subplot(1,1,1)
ax.plot(j, kern, 'b-', label=r'$\mathrm{kernel}$')
ax.plot(j, kern_inv, 'r-', label=r'$\mathrm{inverse \ kernel}$')
ax.legend(loc='upper right', frameon=False)
ax.set_xlabel(r'$\mathrm{displacement}$', fontsize=14)

# Generate covariance matrix from kernel
cov = np.zeros((j.size, j.size))
cov_inv = np.zeros((j.size, j.size))

for i in xrange(j.size):
    cov[i,:] = np.roll(kern, i+j.size/2+1)
    cov_inv[i,:] = np.roll(kern_inv, i+j.size/2+1)

# Verify that the we've calculated the correct inverse
I = np.dot(cov, cov_inv)

# Plot covariance matrix, its inverse and their matrix product
fig = plt.figure(figsize=(12,4), dpi=150)

ax = fig.add_subplot(1,3,1)
ax.imshow(cov.T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='binary', vmin=0., vmax=1.)
ax.set_title(r'$\Sigma$', fontsize=14)
ax.set_ylim(ax.get_ylim()[::-1])

ax = fig.add_subplot(1,3,2)
gamma = 2.
img = np.power(np.abs(cov_inv), 1./gamma) * np.sign(cov_inv)
vmax = np.max(np.abs(img))
ax.imshow(img.T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='bwr_r', vmin=-vmax, vmax=vmax)
ax.set_title(r'$\Sigma^{-1}$', fontsize=14)
ax.set_ylim(ax.get_ylim()[::-1])

ax = fig.add_subplot(1,3,3)
gamma = 2.
img = np.power(np.abs(I), 1./gamma) * np.sign(I)
vmax = np.max(np.abs(img))
ax.imshow(img.T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='bwr_r', vmin=-vmax, vmax=vmax)
ax.set_title(r'$\Sigma^{-1} \Sigma$', fontsize=14)
ax.set_ylim(ax.get_ylim()[::-1])

plt.show()

cov_inv_2D.py

import numpy as np
import matplotlib.pyplot as plt

# Properties of kernel
a = 0.     # Coeff. on Gaussian term (exponential term gets 1-a)
ell = 1.    # Exponential correlation length
sigma = 1.  # Gaussian std. dev. parameter

# 1D pixel array
i = np.arange(-30,31,1)
j = np.arange(-30,31,1)

ii,jj = np.meshgrid(i,j)

# Kronecker Delta
d = np.zeros(ii.shape)
d[-1,-1] = 1.

# Kernel
dist = np.sqrt((ii)**2. + jj**2.)
kern_gauss = np.exp(-0.5 * (np.abs(dist)/sigma)**2.)
kern_exp = np.exp(-np.abs(dist/ell))
kern = a*kern_gauss + (1.-a)*kern_exp

# Inverse Kernel from Fourier-space formula: K^-1 = d/K
kern_inv = np.real(np.fft.ifft2(np.fft.fft2(d) / np.fft.fft2(kern)))

# Plot kernel and its inverse
fig = plt.figure(figsize=(8,4), dpi=150)

i0, i1 = i.size/2-15, i.size/2+16   # Only show the central part of the kernel

ax = fig.add_subplot(1,2,1)
ax.imshow(kern[i0:i1,i0:i1].T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='binary', vmin=0., vmax=1.)

ax = fig.add_subplot(1,2,2)
gamma = 1.
img = np.power(np.abs(kern_inv), 1./gamma) * np.sign(kern_inv)
vmax = np.max(np.abs(img))
ax.imshow(img[i0:i1,i0:i1].T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='bwr_r', vmin=-vmax, vmax=vmax)

# Generate covariance matrix from kernel
n_pix = i.size * j.size
cov = np.zeros((n_pix, n_pix))
cov_inv = np.zeros((n_pix, n_pix))

kern = kern.flatten()
kern_inv = kern_inv.flatten()

for l in xrange(n_pix):
    cov[l,:] = np.roll(kern, l+n_pix/2+1)
    cov_inv[l,:] = np.roll(kern_inv, l+n_pix/2+1)

# Verify that the we've calculated the correct inverse
I = np.dot(cov, cov_inv)

# Plot covariance matrix, its inverse and their matrix product
fig = plt.figure(figsize=(12,4), dpi=150)

k0 = n_pix/2 - 50
k1 = n_pix/2 + 50

ax = fig.add_subplot(1,4,1)
ax.imshow(cov.T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='binary', vmin=0., vmax=1.)
ax.set_title(r'$\Sigma$', fontsize=14)
ax.set_ylim(ax.get_ylim()[::-1])

ax = fig.add_subplot(1,4,2)
gamma = 2.
img = np.power(np.abs(cov_inv), 1./gamma) * np.sign(cov_inv)
vmax = np.max(np.abs(img))
ax.imshow(img.T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='bwr_r', vmin=-vmax, vmax=vmax)
ax.set_title(r'$\Sigma^{-1}$', fontsize=14)
ax.set_ylim(ax.get_ylim()[::-1])

ax = fig.add_subplot(1,4,3)
gamma = 2.
img = np.power(np.abs(I), 1./gamma) * np.sign(I)
vmax = np.max(np.abs(img))
ax.imshow(img.T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='bwr_r', vmin=-vmax, vmax=vmax)
ax.set_title(r'$\Sigma^{-1} \Sigma$', fontsize=14)
ax.set_ylim(ax.get_ylim()[::-1])

ax = fig.add_subplot(1,4,4)
gamma = 2.
img = np.power(np.abs(I), 1./gamma) * np.sign(I)
vmax = np.max(np.abs(img))
ax.imshow(img[k0:k1,k0:k1].T, origin='lower', aspect='equal',
          interpolation='nearest', cmap='bwr_r', vmin=-vmax, vmax=vmax)
ax.set_title(r'$\Sigma^{-1} \Sigma \ \left( \mathrm{zoom \ in} \right)$', fontsize=14)
ax.set_ylim(ax.get_ylim()[::-1])

np.fill_diagonal(I, 0.)

pctiles = [0., 0.5, 1., 10., 25., 50., 75., 90., 95., 99., 99.5, 99.9, 99.99, 100.]
v_pctiles = np.percentile(I, pctiles)

print 'Percentiles of off-diagonal terms in cov^-1 cov'
print '==============================================='

for p,v in zip(pctiles, v_pctiles):
    print r'{:.2f}% : {:.4f}'.format(p,v)

plt.show()