# Chapter 2 
> Linear Algebra and Machine Learning

## Lecture 9
___

### Review of Linear Algebra

Reference Books: Matrix Cookbook by Kaare Brandt Petersen & Michael Syskind Pedersen, 2012



$A \in \mathbb{R}^{n \times m}, n\text{ rows and } m\text{ columns}$ 

range($A$):=span{\underline{a}$_1$,...,\underline{a}$_m$}

null($A$):={\underline{x} $\in \mathbb{R}^{m}$|$A$\underline{x}=0}

Column rank = Row rank = number of linearly independent vectors.

Full rank: rank(A) = min{m,n}.

A \textbf{nonsingular} or \textbf{ invertible} matrix is a square matrix of full rank

Angle between two vectors: $\displaystyle \alpha = cos^{-1}(\frac{\mathbf{x}^T\mathbf{y}}{||\mathbf{x}||\cdot ||\mathbf{y}|| })$

$Q$ matrix is \textbf{ unitary} or \textbf{ orthogonal} if $Q^T=Q^{-1}$:
$||Q\mathbf{x}||=||\mathbf{x}||$ rotation

\textbf{ Norms}:

$\displaystyle ||\mathbf{x}||_p:=(\sum_{j=1}^n|\mathbf{x}_j|^p)^{1/p}$

$\displaystyle ||\mathbf{x}||_1:=\sum_{j=1}^n|\mathbf{x}_j|$

$\displaystyle ||\mathbf{x}||_2:=(\mathbf{x}^T\mathbf{x})^{1/2}$

$\displaystyle ||\mathbf{x}||_\infty:=\max_{1\le j\le n}|\mathbf{x}_j|$

$\displaystyle 
||A||_{(m,n)}:=\sup_{x\in \mathbf{R}^m, \mathbf{x} \ne 0}\frac{||A\mathbf{x}||_n}{||\mathbf{x}||_m}=\sup_{||\mathbf{x}||_m=1}||A\mathbf{x}||_n$

$\displaystyle ||A||_1 = \text{ max column in } A$ <br>
$\displaystyle ||A||_\infty = \text{ max row in } A$

Frobenins (Hilbert-Schmidt norm): 

$\displaystyle ||A||_F = (\sum_{i=1}^n\sum_{j=1}^m |a_{i,j}|^2)^{1/2}=\sqrt{\text{ Tr}(A^TA)}=\sqrt{\text{Tr}(AA^T)}\text{ , where Tr}(B) = \sum_j b_{jj}, \text{ sum of diagonal entries.}$

$||QA||_2=||A||_2 \,\&\, ||QA||_F=||A||_F $

\textbf{Singular Value Decomposition} (SVD): rotation & stretching of a basis

$U$ - left singular matrix, $V$- right singular matrix, $\Sigma$ is the diagonal entry matrix.



\begin{align}
\text{To find the } V, \text{we note that } A^TA &= V\Sigma U^{-1}U\Sigma V^T \\
& = V\Sigma^2V^T \\
&= V\left[ 
\begin{array}{ccc}
\sigma_1^2 &  & \\
& \sigma_2^2 &  \\
& & ... \\
& & \sigma_n^2  \\
\end{array}
\right]V^T
\end{align}

The eigenvectors of this matrix $A^TA$ will give us the vectors $\mathbf{v}_i$, and the eigenvalues will give the numbers $\sigma_i$. 

Similarly, the matrix $AA^T$ gives us info for $U$. Examples can be seen here [SVD examples](
https://ocw.mit.edu/courses/mathematics/18-06sc-linear-algebra-fall-2011/positive-definite-matrices-and-applications/singular-value-decomposition/MIT18_06SCF11_Ses3.5sum.pdf#:~:text=The%20singular%20value%20decomposition%20of%20a%20matrix%20is,its%20eigenvectors%20are%20orthogonal%20and%20we%20can%20write)

<div>

</div>




## Lecture 10
___

### PCA (principal component analysis)

$\overline{x} $ center data by subtraction the mean

$\frac{1}{N-1}XX^T$ is simply the empirical approximation to Cov(X) which henceforth we denote as $C_x$. $C_x$ is non-negative definite & Symmetric (NDS). Hence it has an eigen decomposition $C_x = Q\Lambda Q^T$. The eigen-vectors of $C_x$ are called the principal components or the Rarhumen-Loere modes or PCA modes of $C_x$.

At the same time write $X=U\Sigma V^T$ then, 
\begin{align}
C_x=\frac{1}{N-1}XX^T&=\frac{1}{N-1}U\Sigma V^TV\Sigma U^T\\
&=\frac{1}{N-1}U\Sigma^2U^T
\end{align}


Thus, columns of U,\textbf{ the left singular vectors of X are precisely the Principal Components of} $C_x$

But why do we care about the PCA?

A random varibale is called Gaussian if its PDF has the form 
$$\Pi(x)=\frac{1}{\sqrt{(2\pi)^ddet(C)}}exp(-\frac{1}{2}(x-m)^TC^{-1}(x-m))$$, where m is mean

\begin{lemma}
Suppose $x \sim \mathcal{N}(m,C),m \in \mathbb{R}^d, C\in \mathbb{R}^{d\times{d}}. \text{ Let } b\in \mathbb{R}^n \& A\in \mathbb{R}^{n\times d} \text{ then } z=Ax+b \text{ is also Gaussian & } z\in \mathcal{N}(b+Am, ACA^T)  $
\end{lemma}

PCA explains the covariance and directions of maximum singular values in our dataset. It is good for low-dimension compression. 


## Lecture 11
___
### From PCA & SVD to Proper Orthogonal Decomposition & Dynamic Mode Decoposition

POD(proper orthogonal decomposition)

$\displaystyle f(x,t):=\sum_{k=0}^{\infty}c_k(t)\psi_k(x)$

In the context, $\psi_k(x)$ are called the POD modes & $c_k(t)$ are called the POD coefficient. Luckily, there is an easy way to compute POD using SVD.

Suppose dynamics of f(x,t) are observed over u discrete set, ie (x_j,t_k) ,j=0...N-1, k=0,...,T-1. Now put it into a matrix D:

$$D:=\begin{bmatrix} 
\\
	d_0=f(x_j,t_0) \;\; | \;\; d_1=f(x_j,t_1) \;\; | \;\; ....\;\; |\;\; d_{T-1}=f(x_j,t_{T-1})\\
    \\
	\end{bmatrix}$$
    
Then, compute the SVD of D.

$D=U\Sigma V^T$, $U$ doesn't change even if the matrix is shuffled around, the direction of time is not important. Then the columns of $U$ are precisely the POD modes & the rows of $\Sigma V^T$ are precisely the POD coefficients.

\textbf{ The power of this method is also important that it could be used to predict the future behavior based on the DMD}. Book:Data-Driven Science and Engineering by Steven




```python
#import numpy as np
#import matplotlib.pyplot as plt
#
# first load frames of cylinder flow simulation
#from google.colab import drive
# drive.mount('/content/drive')
#
data = np.load(
    "/content/drive/MyDrive/Courses/AMATH482582-WIN2022/Notebooks/CylData/cyldata.npy")
Nx = 200
Ny = 50
```


```python
# ax0 is the flattened frames and ax1 are the time frames

# compute svd of data

import numpy.matlib as matlib

mean_data = np.mean(data, 1)

centered_data = data - np.transpose(matlib.repmat(mean_data, 200, 1))

dU, ds, dVt = np.linalg.svd(centered_data, full_matrices=False)
```


```python
# plot some of the singular values and some of the principal modes


fig, ax = plt.subplots(1, 1, figsize=(8, 8))

ax.plot(np.log(ds[0:100]))
ax.set_xlabel('index $j$')
ax.set_ylabel('$\log(\sigma_j)$')
#ax.set_xlim(0, 150)


fig, ax = plt.subplots(8, 2, figsize=(20, 20))

for j in range(8):

    ax[j][0].imshow(np.reshape(dU[:, j], (Ny, Nx)), cmap='bwr')
    ax[j][0].get_xaxis().set_visible(False)
    ax[j][0].get_yaxis().set_visible(False)
    ax[j][0].set_aspect('equal')

    ax[j][1].plot(np.abs(dVt[j, :]))
    ax[j][1].set_xlabel('Time step')

    if j == 0:
        ax[j][0].set_title('POD Modes', fontsize=30)
        ax[j][1].set_title('|POD Coeff.|', fontsize=30)

plt.show()
```


```python
# approximate the dynamics using the POD modes only keeping the first 40 and compare to original data.

dss = np.copy(ds)

dss[10:None] = 0

# low rank approx
approx_centered_data = np.dot(dU, np.dot(np.diag(dss), dVt))

print(approx_centered_data.shape)

# add the mean back
approx_data = approx_centered_data + \
    np.transpose(matlib.repmat(mean_data, 200, 1))
```


```python
# side by side comparison of original and approximate dynamics

fig, ax = plt.subplots(4, 2, figsize=(20, 15))

frm_indx = [10, 100, 150, 199]

for j in range(4):

    ax[j][0].imshow(np.reshape(data[:, frm_indx[j]], (Ny, Nx)), cmap='bwr')
    ax[j][0].get_xaxis().set_visible(False)
    ax[j][0].get_yaxis().set_visible(False)
    ax[j][0].set_aspect('equal')

    ax[j][1].imshow(np.reshape(
        approx_data[:, frm_indx[j]], (Ny, Nx)), cmap='bwr')
    ax[j][1].get_xaxis().set_visible(False)
    ax[j][1].get_yaxis().set_visible(False)
    ax[j][1].set_aspect('equal')

    if j == 0:
        ax[j][0].set_title('Original', fontsize=30)
        ax[j][1].set_title('POD Approximation', fontsize=30)
```

### DMD (Dynamic Modes Decomposition) - developed 2010s, some work done at UW:

Look back to the POD, the POD modes are smooth & structured, the POD coefficients are rough & chaotic. This is because POD only sees spatial structure & ignores temporal dynamics. DMD attempts to address this issue by modeling the dynamics of the data:

$$d_{k+1}=Ad_k$$
$$f(x,t_{k+1}){\leftarrow}A\leftarrow f(x,t_k)$$

Lets reconstruct our data as follows, one keeping the T-1 data from 0 to T-2, the other keeping the T-1 data from 1 to T-1.

$$D_1:=\begin{bmatrix} 
\\
	d_0 \;\; | \;\; d_1\;\; | \;\; ....\;\; |\;\; d_{T-2}\\
    \\
	\end{bmatrix}$$

$$D_2:=\begin{bmatrix} 
\\
	d_1 \;\; | \;\; d_2\;\; | \;\; ....\;\; |\;\; d_{T-1}\\
    \\
	\end{bmatrix}$$
    
Then DMD seeks to find a matrix $A$ such that 
$$D_2 \approx AD_1 $$

This approx is then done in a best-fit 


$$A=\min_{B\in \mathbb{R}^{N\times N},\text{ rank(B)}\le r}||BD_1-D_2||_F^2$$

So what does DMD actually do?

\begin{align}
d_k&=Ad_{k-1}=A(Ad_{k-2})=...\\
d_k&=A^{k-1}d_0=Q\Lambda^{k-1}Q^{-1}d_0 \\
&=Q\Lambda^{k-1}b=\sum_{j=1}^rq_{j}\lambda_{j}^{k-1}b_j
\end{align}


```python
from pydmd import DMD
```


```python
# first we create a DMD object
dmd = DMD(svd_rank=10)  # a rank 10 approximation to the dynamic matrix A

# simply pass centered data set of snapshots to the dmd class
dmd.fit(centered_data)

dmd.plot_eigs(show_axes=True)
```


```python
# now plot the DMD modes and the temporal dynamics as we did for POD

fig, ax = plt.subplots(8, 4, figsize=(30, 18))

for j in range(8):

    ax[j][0].imshow(np.reshape(np.real(dmd.modes[:, j]), (Ny, Nx)), cmap='bwr')
    ax[j][0].get_xaxis().set_visible(False)
    ax[j][0].get_yaxis().set_visible(False)
    ax[j][0].set_aspect('equal')

    ax[j][1].plot(np.real(dmd.dynamics[j, :]))
    ax[j][1].set_xlabel('Time step')

    ax[j][2].imshow(np.reshape(np.imag(dmd.modes[:, j]), (Ny, Nx)), cmap='bwr')
    ax[j][2].get_xaxis().set_visible(False)
    ax[j][2].get_yaxis().set_visible(False)
    ax[j][2].set_aspect('equal')

    ax[j][3].plot(np.imag(dmd.dynamics[j, :]))
    ax[j][3].set_xlabel('Time step')

    if j == 0:
        ax[j][0].set_title('Real DMD Modes', fontsize=30)
        ax[j][1].set_title('Real DMD Dynamics', fontsize=30)
        ax[j][2].set_title('Imag DMD Modes', fontsize=30)
        ax[j][3].set_title('Imag DMD Dynamics', fontsize=30)


plt.show()
```

 ## Lecture 12
 ___
 ### Introduction to Machine Learning
 
 Book: The Elements of Statistical Learning by Hastie
 
 

\textbf{Supervised Learning}: With labels.
Predict/classify data given a training dataset. Eg: REgression, classification, function approximation, etc.

\textbf{Un-supervised Learning}:
No labels.
Find meaningful structure in dataset. Eg: Clustering, dimensionality reduction, Generation modeling.

### Supervised Learning

The function model assumes there exists a function
$\displaystyle f^+:X\rightarrow y$ so that $y_j=f^+(x_j)+\epsilon_j$, where $\epsilon_j$ are some noise that may be in the output or our observation of the $f^+(x_j)$

By far the most common assumption is Gaussian noise 

$$\epsilon_j \approx \mathcal{N}(0,\sigma^2)$$

This implies $\displaystyle y_i|x_j \approx \mathcal{N}(f^+(x_j),\sigma^2)$ 

$$\Pi(y_j|x_j)\propto exp(-\frac{1}{2\sigma^2}|f^+(x_j)-y_j|^2)$$, $\Pi$ is the PDF of $y$ for fixed $x_j$

For Euclidern norm, this is called a maximum likelihood estimate (MLE):

$$f_{MLE}=arg\min_f \frac{1}{2\sigma^2}||f(X)-Y||^2$$

At this moment, it is useless without a model, since there are many solutions. One of the most simple model is \textbf{ linear regression}.

$$f_{MLE}\equiv\beta_{MLE}=argmin \frac{1}{2\sigma^2}||A\beta-Y||^2$$, where A=$
\left[
\begin{array}{ccc}
1& & x_0^T \\
.& & . \\
.& & . \\
1& & x_{N-1}^T
\end{array}
\right]
$.

Therefore, MLE is nothing but a least square solution to the problem. Typically, the system is over-determined. Solution is given by solving the normal equations,

$$\frac{\partial}{\partial\beta}(\frac{1}{2\sigma^2}(A\beta-y)^T(A\beta-y))=\frac{1}{\sigma^2}A^T(A\beta-y)=0$$

\begin{align}\implies A^T(A\beta-y)&=0\\
\beta&=(A^TA)^{-1}A^Ty
\end{align}



## Lecture 13
___
### Evaluating SL models

\textbf{ Regularization/penalization/shrinkage}:
we consider $\displaystyle \hat{\beta}=argmin \frac{1}{2\sigma^2}||A\beta-y||^2+\frac{\lambda}{2}||\beta||_p^p$

$\lambda\ge 0$ is called the regularization/penalty parameter & $p \ge 1$ denotes the choice. $p=2 $ for Ridge regression. 

$$\beta=(\frac{1}{\sigma^2}A^TA+\lambda I)^{-1}A^Ty$$

So doing SVD of A,

$$\frac{1}{\sigma^2}A^TA+\lambda I=V(\frac{1}{\sigma^2}\Sigma^2+ \lambda I)V^T$$, the diagonals are non-negative, eliminating zeros so A matrix can be invertible.

Again, the choice of $\lambda$ is important for stability and accuracy.

### Training & Testing Errors

\textbf{Training mean squared error (MSE)}

${X,Y}$ - training set, used to find $\hat{f}(\equiv \hat{\beta})$

${X',Y'}$ - testing set, used for validation.

Analyzing the MSE doesn't mean which model is always better. It is still important for choice of $\lambda$


## Lecture 14
___

### Model tuning with Cross Validation

$\lambda$ too small, the model is basically memorizing all the train data. The test error is large. (over-fitting, high-variance)

$\lambda $ too large, model is too simply biased.()

We want the best test error which is the smallest. However, in real life, we don't know the test error.



\textbf{Cross Validation}: Split the train data to k parts, or k-fold.

Randomly permute the data pairs-- $\mathbf{x}=\{x_{10},x_{-1},...,x{_13}\}$ and responding $\mathbf{y}$. 

Then split the data $\mathbf{x} \,\&\, \mathbf{y} $ into K-subsets. Iterate over $k=0,...k=K-1$ and fit the model to the training data with the k-th fold removed. 

Finally, calculate the CV prediction error (CV cost) with changing $\lambda$

$$CV(\hat{f},\lambda) := \frac{1}{N}\sum_{k=0}^{K-1}||\hat{f}(\mathbf{x}_k,\lambda)-\mathbf{y}_k||^2$$



## Lecture 15
___
### Introduction to Kernel Methods

\textbf{Kernel NDS}: We say K is non-negative definite & symmetric (NDS) if 

\begin{itemize}
\item $K(\mathbf{x},\mathbf{x'})=K(\mathbf{x'},\mathbf{x})$
\item For any set of points $(x_0,...x_n)$, the matrix 
$(K)_{ij}=K(x_i,x_j) $ is NDS.
\item $\displaystyle K(\mathbf{x},\mathbf{x'})=\sum_{j=0}^\infty\lambda_j\psi_j(\mathbf{x})\psi_j(\mathbf{x'}) \text{, the number }\lambda_j\ge 0 $ are the eigenvalues, while the $\psi_j$ are eigenfunctions.

\end{itemize}

Mercer Theorem says if a matrix $A$ is NDS then $A=Q\Lambda Q^{-1}$

$$\implies A=\sum_j\lambda_jq_jq_j^T$$

Define the functions,

$$F_j(x)=\sqrt{\lambda_j}\psi_j(x)$$

along with the mapping,

$$F(\mathbf{x})=(F_0(\mathbf{x}),F_1(\mathbf{x}),....)$$

Functions in RKHS(Reproducing Kernel to Hilbert Space) have the nice properties:

$\displaystyle K (\mathbf{x},\mathbf{x'})=\sum_j F_j(\mathbf{x})F_j(\mathbf{x'}) \,\&\, \text{ so, for a fixed }\mathbf{x}$ we have that $\displaystyle f(\mathbf{x})=\sum_j c_jF_j(\mathbf{x})=(\{c_j\}_{j=0}^\infty,\{F_j(\mathbf{x})\}_{j=0}^\infty)_{l^2}=<f,K(\mathbf{x},)>_{H_k}$

### Kernel Interpolation
Suppose we have an interpolation problem, given x datapoints and y datapoints. we wish to find $\displaystyle f=\sum_{j=0}^\infty c_jF_j$ so we need to solve $\displaystyle \sum_{j=0}^\infty c_jF_j（x_i)=y_i, i=0,...,N-1$ We want the solution with minimized $H_k$ norm. Therefore, we wish to have 

$$minimize \sum_{j=0}^\infty c_j^2$$

$$s.t \sum_{j=0}^\infty c_jF_j(x_i)=y_i$$

What it reduces to is a solution of the form 

$$f(\mathbf{x})=\sum_{j=0}^{N-1}a_jK(\mathbf{x_j},\mathbf{x})$$

The interpolation constraints then tell us that $\displaystyle \sum_{j=0}^{N-1}a_jK(\mathbf{x}_j,\mathbf{x}_i)=y_i$

\begin{align}
\quad\quad\quad\quad\quad\quad\quad\quad \mathbf{a}&=(a_0,....,a_{N-1})\\
\implies \Theta \mathbf{a}=\mathbf{y},\quad\quad\mathbf{y}&=(y_0,....y_{N-1}) \\\quad\quad \Theta_{ji}&=K(\mathbf{x}_j,\mathbf{x}_i)
\end{align}

Thus, $\mathbf{a}=\Theta^{-1}\mathbf{y}$, now the matrix is invertible provided that it is NDS, and the $x_j$ are distinct.


```python
import numpy as np
import matplotlib.pyplot as plt
```


```python
# define function to be interpolated 

def f(x):

  val = np.exp( - x**2/0.05 )*np.cos(np.pi*20*x) + 0.5*np.tanh(5*(x - 0.5))

  return val
```


```python
N = 100
N_p = 5

x = np.linspace(0,1, N) # grid used for plotting  

x_p = np.linspace(0,1, N_p) # set of points for interpolation

# plot the function 

fig,ax = plt.subplots(1,1, figsize=(10,8))

ax.plot(x, f(x))
ax.scatter(x_p, f(x_p), color='r')
ax.set_xlabel('x', fontsize=20)
ax.set_ylabel('f(x)', fontsize=20)
```


```python
# define the kernel we wish to use 
# Gaussian Kernel

def K(x1, x2, l):

  val = np.exp( - (np.abs(x1 - x2)**2)/(l**2) )

  return val


# plot K(x, 0.5) for illustration 

l = [0.05, 0.1, 0.25]

fig,ax = plt.subplots(1,1, figsize=(10,8))

for i in range(3):
  
  ax.plot(x, K(x, 0.5, l[i]), label='l = '+str(l[i]))

ax.set_xlabel('x', fontsize=20)
ax.set_ylabel('K(x, 0.5)', fontsize=20)
plt.legend(fontsize=20)
```


```python
# compute kernel interpolands and plot them for each choice of l

# construct y vector, just for uniform notation with notes
y = f(x_p)

Theta = np.zeros( (N_p, N_p) )

fig,ax = plt.subplots(1,1, figsize=(10,8))

ax.plot(x, f(x), label='f')

for i in range(3):
  # construct kernel matrix (maybe not the most efficient way but ok)
  for j in range(N_p):
    Theta[j,:] = K( x_p[j], x_p, l[i] )

  # compute kernel interpolant 

  a = np.linalg.solve(Theta, y)

  # plot the interpolant along with original function and interpolation data

  KV = np.zeros( (N_p, N) )

  for j in range(N_p):
    KV[j, :] = K( x_p[j], x, l[i] )

  f_interp = np.dot(a, KV)

  ax.plot(x, f_interp, label='l = ' + str(l[i]))


ax.scatter(x_p, y, color='r' )

ax.set_xlabel('x', fontsize=20)
plt.legend(fontsize=20)




```

## Lecture 16
___
### Kernel Ridge Regression 

Mercers Theorem:
$$K(\mathbf{x},\mathbf{x'})=\sum_{j=0}^\infty F_j(\mathbf{x})F_j(\mathbf{x'})$$
,where $F_j$ are the features of $K$, $F_j\rightarrow 0 \text{ as } j\rightarrow \infty$. This regression is capable of solving cases with infinite number of features but at the cost of constructing a finite $n\times n$ matrix. 

\textbf{Ridge Regression}: 
$$\min_{\beta}||A\beta-y||^2+\lambda||\beta||^2$$

where A is our feature matrix often of the form

$$A=\left[\begin{array}{ccc}
F_0(x_0) &F_1(x_1)&...\\
F_0(x_1) & F_1(x_1)&...\\
.. & .. \\
F_0(x_{N-1}) &F_1(x_{N-1})&..
\end{array}\right]$$

while we see $||f||^2_{H_k}=\sum_{j=0}^{J-1}\beta_j^2=||\beta||^2$

\textbf{Kernel Ridge Regression}:
$$\min_{f}||f(X)-Y||^2+\lambda||f||_{H_k}^2$$

\textbf{Representer Theorem}:
$$\hat{f}(\mathbf{x})=\sum_{n=0}^{N-1}\hat{a}_nK(\mathbf{x}_n,\mathbf{x})$$

,with $\hat{a}=(\hat{a}_0,...,\hat{a}_{N-1})$ being minimizer of 
$$minimize \;||\Theta\mathbf{a}-Y||^2+\lambda\mathbf{a}^T\Theta\mathbf{a}$$

The solution is found by differentiating it, and it is exactly $\displaystyle \mathbf{a}=(\Theta+\lambda I)^{-1}y$

### Example: House Pricing in Taiwan

Steps to do:

\begin{itemize}
\item Create corresponding train and test data.
\item Normalize the train and test data: (data-mean)/std
\item Choose and fit Kernel functions (example: "rbf"-Gaussian Kernel)
\item Create discrete $\lambda \text{ and } \gamma $ parameters. 
\item Double for-loop to calculate cross-validation scores (scoring="neg_mean_squared_error") and store the scores. 
\item Find the best score for choosing parameters
\item Predict
\end{itemize}


## Lecture 17
___
### Kernel Perspective vs. Feature Perspective

FR: $minimize \;||A\beta-Y||^2+\lambda||\beta||^2$ \\

KR: $minimize \; ||\Theta a-Y||^2 +\lambda a^T\Theta a$

where $A_{ij}=F_j(x_i)$ matrix of feature values,
$\Theta_{ij}=K(x_i.x_j)$

The size of computation is different in that, provided a large data set with data points far more than feature points, we might want to use FR because KR would generate a matrix of NXN dimension which costs too much.

For KR, we have infinite many features and more freedom to choose.

### Connecting Kernels to Fourier Series & Wavelets

From discrete Fourier Series we note that $\displaystyle f(t)=\sum_{k=-N/2}^{N/2-1}\hat{f}_kexp(ikt)$. We than identify a kernel called the \textbf{ Dirichlet Kernel}:
$$K(t,s)=\sum_{k=-N/2}^{N/2-1}exp(ik(t-s))$$
This gives the identity that 

$$\int_0^{2\pi}K(t,s)f(s)\,ds=\sum_{k=-N/2}^{N/2-1}\hat{f}_kexp(ikt)$$

Then, we also note that linear projection is not the only transformation, non-linear is also possible.



8266bbc8300edcc12b57efbf52ed7806fc99b3b3
ipynb
Jupyter Notebook
course_notes/.ipynb_checkpoints/Chapter2-checkpoint.ipynb
raph651/Amath-582-Data-Analysis
c1d72d897b7611652c7fe1f71c5439062b8bdf9e
MIT
course_notes/.ipynb_checkpoints/Chapter2-checkpoint.ipynb
raph651/Amath-582-Data-Analysis
c1d72d897b7611652c7fe1f71c5439062b8bdf9e
MIT
course_notes/.ipynb_checkpoints/Chapter2-checkpoint.ipynb
raph651/Amath-582-Data-Analysis
c1d72d897b7611652c7fe1f71c5439062b8bdf9e
MIT
Qwen/Qwen-72B
1. YES
2. YES