In [9]:
import numpy as np
import scipy.special as ss
import matplotlib.pyplot as plt
In [10]:
def kernel(x, N=5, spectrum=None):
    """
    Evaluate the reproducing kernel.

    Parameters
    ----------
    x : array_like
        Points at which to evaluate the kernel.
    N : int, optional
        Kernel order.  Defaults to 5.  If the spectrum
        is specified, the kernel order is not used.
    spectrum : array_like, optional
        Weights for each polynomial component of the kernel.
        Ones by default.

    """
    out = np.zeros_like(x)
    
    if spectrum is None:
        spectrum = np.ones(N + 1)
        
    N = len(spectrum)
    
    for l in range(N):
        Pl = ss.legendre(l)
        out += spectrum[l] * (2 * l + 1) / (4 * np.pi) * Pl(x)
    return out
In [11]:
u = np.linspace(-1, 1, 1000)
spectrum = np.ones(5)
k = kernel(u, spectrum=spectrum)

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,4))

ax1.plot(u, k)
ax1.set_title('Kernel value over the interval [-1, 1]')

ax2.plot(spectrum, 'o-')
ax2.set_title('Extended kernel spectrum')
ax2.set_ylim(-0.1, 1.1);
In [12]:
sqrt(2)
Out[12]:
1.4142135623730951
In [13]:
theta = np.linspace(0, 2 * np.pi, 1000)
ang = np.cos(theta)
n1, n2 = 2, 3
k1 = kernel(ang, n1)
k2 = kernel(ang, n2)
f, ax = subplots(subplot_kw=dict(polar=True))
ax.plot(theta, abs(k1), color='red',label=str(n1))
ax.plot(theta, abs(k2), label=str(n2))
ax.legend()
ax.set_title('Kernel magnitude with varying angle');
In [14]:
spectrum = np.zeros(10)
spectrum[:5] = 1
spectrum[5:] = np.linspace(1, 0, 5)
k = kernel(u, spectrum=spectrum)

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,4))

ax1.plot(u, k)
ax1.set_title('Linearly extended kernel')

ax2.plot(spectrum, 'o-')
ax2.set_title('Extended kernel spectrum')
ax2.set_ylim(-0.1, 1.1);
In [15]:
import scipy.optimize as opt

u0 = 0.8
u = np.linspace(-1 + 1e-10, u0, 200)

def full_spectrum(partial_spectrum):
    return np.hstack((np.ones(5), partial_spectrum))

def cost(spectrum):
    return np.sum(np.abs(kernel(u, spectrum=full_spectrum(spectrum))))

opt_spectrum = full_spectrum(opt.fmin(cost, np.ones(5)))
print opt_spectrum

f, (ax1, ax2) = plt.subplots(1, 2)

u = np.linspace(-1, 1, 1000)
ax1.plot(u, kernel(u, spectrum=opt_spectrum))
ax1.set_title('Optimally extended kernel')

ax2.plot(opt_spectrum, 'o-')
ax2.set_title('Kernel spectrum')
ax2.set_ylim(-0.1, 1.1);
Optimization terminated successfully.
         Current function value: 6.628413
         Iterations: 204
         Function evaluations: 327
[ 1.          1.          1.          1.          1.          0.81402883
  0.63712395  0.43044675  0.24297284  0.11031011]
In [15]: