Fixed type and shifted Ising model to end for Regression slides
This commit is contained in:
mhjensen
@@ -1130,7 +1130,7 @@ invert the matrix $\bm{X}^T\bm{X}$. Using a polynomial expansion as we
|
||||
did both for the masses and the fitting of the equation of state,
|
||||
leads to row vectors of the design matrix which are essentially
|
||||
orthogonal due to the polynomial character of our model. Obtaining the inverse of the design matrix is then often done via a so-called LU, QR or Cholesky decomposition.
|
||||
More material to come here.
|
||||
|
||||
|
||||
|
||||
This may
|
||||
@@ -1147,218 +1147,6 @@ inversion algorithm. Thereafter we dive into the math of the SVD.
|
||||
!eblock
|
||||
|
||||
|
||||
# todo: change model here.
|
||||
|
||||
!split
|
||||
===== The Ising model =====
|
||||
|
||||
The one-dimensional Ising model with nearest neighbor interaction, no
|
||||
external field and a constant coupling constant $J$ is given by
|
||||
|
||||
!bt
|
||||
\begin{align}
|
||||
H = -J \sum_{k}^L s_k s_{k + 1},
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
where $s_i \in \{-1, 1\}$ and $s_{N + 1} = s_1$. The number of spins
|
||||
in the system is determined by $L$. For the one-dimensional system
|
||||
there is no phase transition.
|
||||
|
||||
We will look at a system of $L = 40$ spins with a coupling constant of
|
||||
$J = 1$. To get enough training data we will generate 10000 states
|
||||
with their respective energies.
|
||||
|
||||
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
||||
import seaborn as sns
|
||||
import scipy.linalg as scl
|
||||
from sklearn.model_selection import train_test_split
|
||||
import tqdm
|
||||
sns.set(color_codes=True)
|
||||
cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
|
||||
|
||||
L = 40
|
||||
n = int(1e4)
|
||||
|
||||
spins = np.random.choice([-1, 1], size=(n, L))
|
||||
J = 1.0
|
||||
|
||||
energies = np.zeros(n)
|
||||
|
||||
for i in range(n):
|
||||
energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))
|
||||
!ec
|
||||
|
||||
Here we use ordinary least squares
|
||||
regression to predict the energy for the nearest neighbor
|
||||
one-dimensional Ising model on a ring, i.e., the endpoints wrap
|
||||
around. We will use linear regression to fit a value for
|
||||
the coupling constant to achieve this.
|
||||
|
||||
!split
|
||||
===== Reformulating the problem to suit regression =====
|
||||
|
||||
A more general form for the one-dimensional Ising model is
|
||||
|
||||
!bt
|
||||
\begin{align}
|
||||
H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
Here we allow for interactions beyond the nearest neighbors and a state dependent
|
||||
coupling constant. This latter expression can be formulated as
|
||||
a matrix-product
|
||||
!bt
|
||||
\begin{align}
|
||||
\bm{H} = \bm{X} J,
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
where $X_{jk} = s_j s_k$ and $J$ is a matrix which consists of the
|
||||
elements $-J_{jk}$. This form of writing the energy fits perfectly
|
||||
with the form utilized in linear regression, that is
|
||||
|
||||
!bt
|
||||
\begin{align}
|
||||
\bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon},
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
We split the data in training and test data as discussed in the previous example
|
||||
|
||||
!bc pycod
|
||||
X = np.zeros((n, L ** 2))
|
||||
for i in range(n):
|
||||
X[i] = np.outer(spins[i], spins[i]).ravel()
|
||||
y = energies
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== Linear regression =====
|
||||
|
||||
In the ordinary least squares method we choose the cost function
|
||||
|
||||
!bt
|
||||
\begin{align}
|
||||
C(\bm{X}, \bm{\beta})= \frac{1}{n}\left\{(\bm{X}\bm{\beta} - \bm{y})^T(\bm{X}\bm{\beta} - \bm{y})\right\}.
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
We then find the extremal point of $C$ by taking the derivative with respect to $\bm{\beta}$ as discussed above.
|
||||
This yields the expression for $\bm{\beta}$ to be
|
||||
|
||||
!bt
|
||||
\[
|
||||
\bm{\beta} = \frac{\bm{X}^T \bm{y}}{\bm{X}^T \bm{X}},
|
||||
\]
|
||||
!et
|
||||
|
||||
which immediately imposes some requirements on $\bm{X}$ as there must exist
|
||||
an inverse of $\bm{X}^T \bm{X}$. If the expression we are modeling contains an
|
||||
intercept, i.e., a constant term, we must make sure that the
|
||||
first column of $\bm{X}$ consists of $1$. We do this here
|
||||
|
||||
!bc pycod
|
||||
X_train_own = np.concatenate(
|
||||
(np.ones(len(X_train))[:, np.newaxis], X_train),
|
||||
axis=1
|
||||
)
|
||||
X_test_own = np.concatenate(
|
||||
(np.ones(len(X_test))[:, np.newaxis], X_test),
|
||||
axis=1
|
||||
)
|
||||
!ec
|
||||
|
||||
!bc pycod
|
||||
def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:
|
||||
return scl.inv(x.T @ x) @ (x.T @ y)
|
||||
beta = ols_inv(X_train_own, y_train)
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Singular Value decomposition =====
|
||||
|
||||
Doing the inversion directly turns out to be a bad idea since the matrix
|
||||
$\bm{X}^T\bm{X}$ is singular. An alternative approach is to use the _singular
|
||||
value decomposition_. Using the definition of the Moore-Penrose
|
||||
pseudoinverse we can write the equation for $\bm{\beta}$ as
|
||||
|
||||
!bt
|
||||
\[
|
||||
\bm{\beta} = \bm{X}^{+}\bm{y},
|
||||
\]
|
||||
!et
|
||||
|
||||
where the pseudoinverse of $\bm{X}$ is given by
|
||||
|
||||
!bt
|
||||
\[
|
||||
\bm{X}^{+} = \frac{\bm{X}^T}{\bm{X}^T\bm{X}}.
|
||||
\]
|
||||
!et
|
||||
|
||||
Using singular value decomposition we can decompose the matrix $\bm{X} = \bm{U}\bm{\Sigma} \bm{V}^T$,
|
||||
where $\bm{U}$ and $\bm{V}$ are orthogonal(unitary) matrices and $\bm{\Sigma}$ contains the singular values (more details below).
|
||||
where $X^{+} = V\Sigma^{+} U^T$. This reduces the equation for
|
||||
$\omega$ to
|
||||
!bt
|
||||
\begin{align}
|
||||
\bm{\beta} = \bm{V}\bm{\Sigma}^{+} \bm{U}^T \bm{y}.
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
Note that solving this equation by actually doing the pseudoinverse
|
||||
(which is what we will do) is not a good idea as this operation scales
|
||||
as $\mathcal{O}(n^3)$, where $n$ is the number of elements in a
|
||||
general matrix. Instead, doing $QR$-factorization and solving the
|
||||
linear system as an equation would reduce this down to
|
||||
$\mathcal{O}(n^2)$ operations.
|
||||
|
||||
|
||||
!bc pycod
|
||||
def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:
|
||||
u, s, v = scl.svd(x)
|
||||
return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y
|
||||
!ec
|
||||
|
||||
!bc pycod
|
||||
beta = ols_svd(X_train_own,y_train)
|
||||
!ec
|
||||
|
||||
When extracting the $J$-matrix we need to make sure that we remove the intercept, as is done here
|
||||
|
||||
!bc pycod
|
||||
J = beta[1:].reshape(L, L)
|
||||
!ec
|
||||
|
||||
A way of looking at the coefficients in $J$ is to plot the matrices as images.
|
||||
|
||||
|
||||
!bc pycod
|
||||
fig = plt.figure(figsize=(20, 14))
|
||||
im = plt.imshow(J, **cmap_args)
|
||||
plt.title("OLS", fontsize=18)
|
||||
plt.xticks(fontsize=18)
|
||||
plt.yticks(fontsize=18)
|
||||
cb = fig.colorbar(im)
|
||||
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
|
||||
plt.show()
|
||||
!ec
|
||||
It is interesting to note that OLS
|
||||
considers both $J_{j, j + 1} = -0.5$ and $J_{j, j - 1} = -0.5$ as
|
||||
valid matrix elements for $J$.
|
||||
In our discussion below on hyperparameters and Ridge and Lasso regression we will see that
|
||||
this problem can be removed, partly and only with Lasso regression.
|
||||
|
||||
In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD?
|
||||
|
||||
|
||||
!split
|
||||
===== Linear Regression Problems =====
|
||||
@@ -1391,7 +1179,7 @@ column vectors. Hence, the rank of $\mathbf{X}$ is equal to the number
|
||||
of linearly independent columns. In this particular case the matrix has rank 2.
|
||||
|
||||
Super-collinearity of an $(n \times p)$-dimensional design matrix $\mathbf{X}$ implies
|
||||
that the inverse of the matrix $\bm{X}^T\bm{x}$ (the matrix we need to invert to solve the linear regression equations) is non-invertible. If we have a square matrix that does not have an inverse, we say this matrix singular. The example here demonstrates this
|
||||
that the inverse of the matrix $\bm{X}^T\bm{X}$ (the matrix we need to invert to solve the linear regression equations) is non-invertible. If we have a square matrix that does not have an inverse, we say this matrix singular. The example here demonstrates this
|
||||
!bt
|
||||
\begin{align*}
|
||||
\bm{X} & = \left[
|
||||
@@ -1527,7 +1315,7 @@ In the general case, where our design matrix $\bm{X}$ has dimension
|
||||
$n\times p$, the matrix is thus decomposed into an $n\times n$
|
||||
orthogonal matrix $\bm{U}$, a $p\times p$ orthogonal matrix $\bm{V}$
|
||||
and a diagonal matrix $\bm{\Sigma}$ with $r=\mathrm{min}(n,p)$
|
||||
singular values $\sigma_i\lg 0$ on the main diagonal and zeros filling
|
||||
singular values $\sigma_i\geq 0$ on the main diagonal and zeros filling
|
||||
the rest of the matrix. There are at most $p$ singular values
|
||||
assuming that $n > p$. In our regression examples for the nuclear
|
||||
masses and the equation of state this is indeed the case, while for
|
||||
@@ -3388,6 +3176,220 @@ plt.show()
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== The Ising model =====
|
||||
|
||||
The one-dimensional Ising model with nearest neighbor interaction, no
|
||||
external field and a constant coupling constant $J$ is given by
|
||||
|
||||
!bt
|
||||
\begin{align}
|
||||
H = -J \sum_{k}^L s_k s_{k + 1},
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
where $s_i \in \{-1, 1\}$ and $s_{N + 1} = s_1$. The number of spins
|
||||
in the system is determined by $L$. For the one-dimensional system
|
||||
there is no phase transition.
|
||||
|
||||
We will look at a system of $L = 40$ spins with a coupling constant of
|
||||
$J = 1$. To get enough training data we will generate 10000 states
|
||||
with their respective energies.
|
||||
|
||||
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
||||
import seaborn as sns
|
||||
import scipy.linalg as scl
|
||||
from sklearn.model_selection import train_test_split
|
||||
import tqdm
|
||||
sns.set(color_codes=True)
|
||||
cmap_args=dict(vmin=-1., vmax=1., cmap='seismic')
|
||||
|
||||
L = 40
|
||||
n = int(1e4)
|
||||
|
||||
spins = np.random.choice([-1, 1], size=(n, L))
|
||||
J = 1.0
|
||||
|
||||
energies = np.zeros(n)
|
||||
|
||||
for i in range(n):
|
||||
energies[i] = - J * np.dot(spins[i], np.roll(spins[i], 1))
|
||||
!ec
|
||||
|
||||
Here we use ordinary least squares
|
||||
regression to predict the energy for the nearest neighbor
|
||||
one-dimensional Ising model on a ring, i.e., the endpoints wrap
|
||||
around. We will use linear regression to fit a value for
|
||||
the coupling constant to achieve this.
|
||||
|
||||
!split
|
||||
===== Reformulating the problem to suit regression =====
|
||||
|
||||
A more general form for the one-dimensional Ising model is
|
||||
|
||||
!bt
|
||||
\begin{align}
|
||||
H = - \sum_j^L \sum_k^L s_j s_k J_{jk}.
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
Here we allow for interactions beyond the nearest neighbors and a state dependent
|
||||
coupling constant. This latter expression can be formulated as
|
||||
a matrix-product
|
||||
!bt
|
||||
\begin{align}
|
||||
\bm{H} = \bm{X} J,
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
where $X_{jk} = s_j s_k$ and $J$ is a matrix which consists of the
|
||||
elements $-J_{jk}$. This form of writing the energy fits perfectly
|
||||
with the form utilized in linear regression, that is
|
||||
|
||||
!bt
|
||||
\begin{align}
|
||||
\bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon},
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
We split the data in training and test data as discussed in the previous example
|
||||
|
||||
!bc pycod
|
||||
X = np.zeros((n, L ** 2))
|
||||
for i in range(n):
|
||||
X[i] = np.outer(spins[i], spins[i]).ravel()
|
||||
y = energies
|
||||
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== Linear regression =====
|
||||
|
||||
In the ordinary least squares method we choose the cost function
|
||||
|
||||
!bt
|
||||
\begin{align}
|
||||
C(\bm{X}, \bm{\beta})= \frac{1}{n}\left\{(\bm{X}\bm{\beta} - \bm{y})^T(\bm{X}\bm{\beta} - \bm{y})\right\}.
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
We then find the extremal point of $C$ by taking the derivative with respect to $\bm{\beta}$ as discussed above.
|
||||
This yields the expression for $\bm{\beta}$ to be
|
||||
|
||||
!bt
|
||||
\[
|
||||
\bm{\beta} = \frac{\bm{X}^T \bm{y}}{\bm{X}^T \bm{X}},
|
||||
\]
|
||||
!et
|
||||
|
||||
which immediately imposes some requirements on $\bm{X}$ as there must exist
|
||||
an inverse of $\bm{X}^T \bm{X}$. If the expression we are modeling contains an
|
||||
intercept, i.e., a constant term, we must make sure that the
|
||||
first column of $\bm{X}$ consists of $1$. We do this here
|
||||
|
||||
!bc pycod
|
||||
X_train_own = np.concatenate(
|
||||
(np.ones(len(X_train))[:, np.newaxis], X_train),
|
||||
axis=1
|
||||
)
|
||||
X_test_own = np.concatenate(
|
||||
(np.ones(len(X_test))[:, np.newaxis], X_test),
|
||||
axis=1
|
||||
)
|
||||
!ec
|
||||
|
||||
!bc pycod
|
||||
def ols_inv(x: np.ndarray, y: np.ndarray) -> np.ndarray:
|
||||
return scl.inv(x.T @ x) @ (x.T @ y)
|
||||
beta = ols_inv(X_train_own, y_train)
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Singular Value decomposition =====
|
||||
|
||||
Doing the inversion directly turns out to be a bad idea since the matrix
|
||||
$\bm{X}^T\bm{X}$ is singular. An alternative approach is to use the _singular
|
||||
value decomposition_. Using the definition of the Moore-Penrose
|
||||
pseudoinverse we can write the equation for $\bm{\beta}$ as
|
||||
|
||||
!bt
|
||||
\[
|
||||
\bm{\beta} = \bm{X}^{+}\bm{y},
|
||||
\]
|
||||
!et
|
||||
|
||||
where the pseudoinverse of $\bm{X}$ is given by
|
||||
|
||||
!bt
|
||||
\[
|
||||
\bm{X}^{+} = \frac{\bm{X}^T}{\bm{X}^T\bm{X}}.
|
||||
\]
|
||||
!et
|
||||
|
||||
Using singular value decomposition we can decompose the matrix $\bm{X} = \bm{U}\bm{\Sigma} \bm{V}^T$,
|
||||
where $\bm{U}$ and $\bm{V}$ are orthogonal(unitary) matrices and $\bm{\Sigma}$ contains the singular values (more details below).
|
||||
where $X^{+} = V\Sigma^{+} U^T$. This reduces the equation for
|
||||
$\omega$ to
|
||||
!bt
|
||||
\begin{align}
|
||||
\bm{\beta} = \bm{V}\bm{\Sigma}^{+} \bm{U}^T \bm{y}.
|
||||
\end{align}
|
||||
!et
|
||||
|
||||
Note that solving this equation by actually doing the pseudoinverse
|
||||
(which is what we will do) is not a good idea as this operation scales
|
||||
as $\mathcal{O}(n^3)$, where $n$ is the number of elements in a
|
||||
general matrix. Instead, doing $QR$-factorization and solving the
|
||||
linear system as an equation would reduce this down to
|
||||
$\mathcal{O}(n^2)$ operations.
|
||||
|
||||
|
||||
!bc pycod
|
||||
def ols_svd(x: np.ndarray, y: np.ndarray) -> np.ndarray:
|
||||
u, s, v = scl.svd(x)
|
||||
return v.T @ scl.pinv(scl.diagsvd(s, u.shape[0], v.shape[0])) @ u.T @ y
|
||||
!ec
|
||||
|
||||
!bc pycod
|
||||
beta = ols_svd(X_train_own,y_train)
|
||||
!ec
|
||||
|
||||
When extracting the $J$-matrix we need to make sure that we remove the intercept, as is done here
|
||||
|
||||
!bc pycod
|
||||
J = beta[1:].reshape(L, L)
|
||||
!ec
|
||||
|
||||
A way of looking at the coefficients in $J$ is to plot the matrices as images.
|
||||
|
||||
|
||||
!bc pycod
|
||||
fig = plt.figure(figsize=(20, 14))
|
||||
im = plt.imshow(J, **cmap_args)
|
||||
plt.title("OLS", fontsize=18)
|
||||
plt.xticks(fontsize=18)
|
||||
plt.yticks(fontsize=18)
|
||||
cb = fig.colorbar(im)
|
||||
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=18)
|
||||
plt.show()
|
||||
!ec
|
||||
It is interesting to note that OLS
|
||||
considers both $J_{j, j + 1} = -0.5$ and $J_{j, j - 1} = -0.5$ as
|
||||
valid matrix elements for $J$.
|
||||
In our discussion below on hyperparameters and Ridge and Lasso regression we will see that
|
||||
this problem can be removed, partly and only with Lasso regression.
|
||||
|
||||
In this case our matrix inversion was actually possible. The obvious question now is what is the mathematics behind the SVD?
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!split
|
||||
===== The one-dimensional Ising model =====
|
||||
|
||||
|
||||
Reference in New Issue
Block a user