packages <- c('tidyverse', 'ggthemes', 'knitr', 'kableExtra', 'Matrix', 'reticulate')
sapply(packages, require, character.only = TRUE, quietly = TRUE)
use_condaenv("r-reticulate")
theme_set(theme_tufte(base_size = 14) +
                theme(panel.border = element_rect(fill = NA),
                      panel.grid.major = element_line(color = "gray78"),
                      legend.background = element_rect(),
                      legend.position = "top",
                      axis.text.x = element_text(angle = 30, hjust = 1),
                      strip.background = element_rect(fill = "grey90", linetype = "blank")))

Problem 1

Here H is always Hilbert space over \(\mathbb R\) with scalar-product \(\langle ., . \rangle\) and \(||.|| = \sqrt{\langle ., . \rangle}\).

1. For a separable Hilbert space H with orthonormal basis \(\mathbf e_n, n \in \mathbb N\), we showed that if \(\lambda \in H^*\) then \(\lambda(\mathbf x) = \langle \mathbf x, \mathbf z \rangle\) for \(\mathbf z = \sum \lambda(\mathbf e_n)\mathbf e_n\). Show that \(\mathbf z \perp \ker (\lambda)\) and that \(||\mathbf z||^2 =\lambda(z)\)

Solution: Since \(\ker (\lambda) = \{\mathbf x \in H, \lambda(\mathbf x) = 0\}\) we can see that \(\lambda(\ker(\lambda)) = \langle \ker(\lambda), \mathbf z \rangle = 0\) which means that \(\mathbf z \perp \ker(\lambda)\).

From the definition of \(\lambda\) we also have that \(\begin{aligned} ||\mathbf z||^2 &= \langle \mathbf z, \mathbf z \rangle = \lambda(z) \end{aligned}\)

2. Show that if \((\mathbf x_k)_{k\in\mathbb N}\) is a Cauchy-sequence in a Hilbert space and \((\mathbf y_k)_{k\in\mathbb N}\) is a subsequence of \((\mathbf x_k)_{k\in\mathbb N}\) i.e. \(\mathbf y_k = \mathbf x_{n_k}\) for all \(k \in \mathbb N\), such that \(lim_{k\rightarrow\infty}\mathbf y_k = x\), then \(lim_{k\rightarrow\infty}\mathbf x_k = x\)

Solution: Let \(\epsilon > 0\). Since \((\mathbf x_k)_{k\in\mathbb N}\) is Cauchy we can take N > 0 such that for all n, m > N we have

\[ ||x_n - x_m|| < \frac\epsilon 2. \]

We know that \((\mathbf y_k)_{k\in\mathbb N}\) is a subsequence and \(lim_{k\rightarrow\infty}\mathbf y_k = x\) so we can take K > 0 suck that for all \(n_k > K\) we have

\[ ||y_k - x|| = ||x_{n_k} - x|| < \frac\epsilon 2. \]

If we now choose M = max(N, K) we have that for all \(n, m, n_k > M\)

\[ ||x_n - x|| \leq ||x_n - x_{n_k}|| + ||x_{n_k} - x|| \leq \frac\epsilon 2 + \frac\epsilon 2 = \epsilon \]

3. Show that if \((\mathbf x_k)_{k\in\mathbb N}\) is a Cauchy-sequence in a Hilbert space, then there is a subsequence \((\mathbf y_k)_{k\in\mathbb N}\) of \((\mathbf x_k)_{k\in\mathbb N}\), such that \(||\mathbf y_{k+1} - \mathbf y_k|| \leq 2^{-k}\) for all \(k \in \mathbb N\)

Solution: Since \((\mathbf x_k)_{k\in\mathbb N}\) is a Cauchy sequence we can create \((\mathbf y_k)_{k\in\mathbb N}\) as a subsequence of \((\mathbf x_k)_{k\in\mathbb N}\) such that it is also Cauchy. Let \(\epsilon > 0\) and choose N > 0 such that \(||y_{k + 1} - y_k|| \leq \epsilon\) for all k > N. In particular, since we can find such N for all \(\epsilon > 0\) we can let \(\epsilon = 2^{-k}\) and thus \(||y_{k + 1} - y_k|| \leq 2^{-k}\)

Problem 2

For a kernel \(K : \mathbb R ^ n \times \mathbb R^n \rightarrow \mathbb R\) and a grid \(X = \{\mathbf x_1, \mathbf x_2, ..., \mathbf x_N\} \subset \mathbb R^n\) we define the matrix \(A_{K,X} := (K(\mathbf x_i, \mathbf x_j))_{i,j=1:N}\)

1. Program the kernel functions \(K_s(\mathbf x, \mathbf y) = \phi_s(||\mathbf x - \mathbf y||)\), where

\[ \begin{aligned} \phi_1(r) &= (1 - r)^4_+(4r + 1) \\ \phi_2(r) &= (1 + r^2)^{-1/2} \\ \phi_3(r) &= \exp(-r^2) \end{aligned} \]

import tensorflow as tf
import numpy as np
def phi_1(x, y):
      input_x = tf.placeholder(dtype = "float32", shape = None)
      input_y = tf.placeholder(dtype = "float32", shape = None)
      r = tf.linalg.norm(input_x - input_y)
      phi = tf.nn.relu(1 - r)**4 * (4 * r + 1)
      with tf.Session() as sess:
            out = sess.run(phi, feed_dict = ({input_x: x, input_y: y}))
      return out
def phi_2(x, y):
      input_x = tf.placeholder(dtype = "float32", shape = None)
      input_y = tf.placeholder(dtype = "float32", shape = None)
      r = tf.linalg.norm(input_x - input_y)
      phi = tf.pow(1 + r**2, -1/2)
      with tf.Session() as sess:
            out = sess.run(phi, feed_dict = ({input_x: x, input_y: y}))
      return out
def phi_3(x, y):
      input_x = tf.placeholder(dtype = "float32", shape = None)
      input_y = tf.placeholder(dtype = "float32", shape = None)
      r = tf.linalg.norm(input_x - input_y)
      phi = tf.exp(-r**2)
      with tf.Session() as sess:
            out = sess.run(phi, feed_dict = ({input_x: x, input_y: y}))
      return out

2. Write a program that given a kernel K and grid \(X = \{\mathbf x_1, \mathbf x_2, ..., \mathbf x_n\} \subset \mathbb R^n\) (vector of vectors), computes the interpolation matrix \(A_{K,X} \in \mathbb R^{N \times N}\)

def interpolate(x, kernel):
      input_x = tf.placeholder("float32", shape = None)
      input_y = tf.placeholder("float32", shape = None)
      r = tf.linalg.norm(input_x - input_y)
      out = np.zeros((x.shape[0], x.shape[0]))
      phi = kernel(r)
      with tf.Session() as sess:
            for i, x1 in enumerate(x):
                  for j, x2 in enumerate(x):
                        out[i, j] = sess.run(phi, feed_dict = ({input_x: x1, input_y: x2}))
      return out
phi_1 = lambda r: tf.nn.relu(1 - r)**4 * (4 * r + 1)
phi_2 = lambda r: tf.pow(1 + r**2, -1./2)
phi_3 = lambda r: tf.exp(-r**2)
x = np.random.normal(0, 1, (10, 10))
print(interpolate(x, phi_1).shape)
## (10, 10)

3. Write a program that given a grid \(X = \{\mathbf x_1, \mathbf x_2, ..., \mathbf x_n\} \subset \mathbb R^n\) (vector of vectors), an interpolation matrix \(A_{K,X} \in \mathbb R^{N \times N}\) and a function \(W:\mathbb R^n \rightarrow \mathbb R\), solves the linear system \(A_{K,X}\mathbf \alpha = \mathbf \beta\), where \(\beta_i = W(\mathbf x_i)\) for \(i = 1: N\)

def kern_solver(x, A, W):
      beta = W(x)
      alpha = np.linalg.solve(A, beta)
      return(alpha)
W = lambda x: np.mean(x, axis = 0)
x = np.random.normal(0, 1, (10, 10))
A = interpolate(x, phi_1)
print(np.linalg.solve(A, W(x)))
## [-0.05590357  0.14193602  0.40256645 -0.2013339  -0.2231713  -0.0501973
##  -0.1598773  -0.07209407  0.46772518 -0.45284432]