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]: