<jupyter_start><jupyter_text>Training a Sequential RBM on a timeseries The purpose of this notebook is to examine how well the sequential Bernoulli RBM fits a temporal series of binary values. The model is explained in a previous[notebook](rbm_models.ipynb), but we shall revisit the mathematics here. For the purposes of the current experiment, we utilise the file `acgl.us.txt` of stock prices, obtained from [Kaggle](https://www.kaggle.com/gulabpatelcovid/stock-market-analysis-and-time-series-prediction/data).<jupyter_code>import sys

sys.path.append("../python")
import os
import numpy as np
import pandas as pd
DATA_PATH = "../../../../projects/financial-timeseries/stocks"  # YMMV
df = pd.read_csv(os.path.join(DATA_PATH, "acgl.us.txt"))
df.head()
df.tail()<jupyter_output><empty_output><jupyter_text>For want of a better idea, we follow the usual idea of tracking the closing price of the stock at the end of each day. Since we currently only have an implementation of a Bernoulli RBM, we shall signal a rise in price (from one day to the next) by 1, anda fall (or the same value) by 0.<jupyter_code>c_series = df.Close.values
d_series = np.diff(c_series)
b_series = np.asarray(d_series > 0, dtype=int)
print(b_series[:10])
print(b_series[-10:])<jupyter_output>[1 0 1 0 1 1 0 0 0 0]
[1 0 0 0 0 1 1 0 1 0]<jupyter_text>Next, we need to form the timeseries into a collection of vectors, based upon sliding a window of fixed length along the series.<jupyter_code>from numpy.lib.stride_tricks import sliding_window_view as _window<jupyter_output><empty_output><jupyter_text>Let us start with the simplest possible Bernoulli RBM, with one input bit and one output bit. We note in advance that the modelled input bit probability, $p(x_1=1)$, must be constant.In fact, the parameter estimates should converge via gradient ascent when the model probability equals the empirical probability (as we shall verify later).<jupyter_code>num_input = 1
num_output = 1
X_train = _window(b_series, num_input)
print(X_train.shape)
print(X_train[:9,:])
print(X_train[-9:,:])
print("There are two possible states:")
idx = [1, 0]
X_states = X_train[idx, 0:1]
print(X_states)
ind0 = X_train[:, 0] == 0
p0 = sum(ind0) / len(ind0)
p1 = sum(~ind0) / len(~ind0)
base_probs = (p0, p1)
print("The empirical probabilities are:")
print("p(x1=0)=%f, p(x1=1)=%f" % base_probs)
from bernoulli_rbm import SequentialBernoulliRBM, ExactSequentialBernoulliRBM
rbm = SequentialBernoulliRBM(
    num_output, num_input, 
)
rbm.fit(X_train)

preds = rbm.reconstruct(X_states)
probs = preds[:, 0]  # p(x_1)
probs = np.array(list(zip(1 - probs, probs)))
print(probs[range(2), X_states[:, 0]])<jupyter_output>[0.47164896 0.52835104]<jupyter_text>Note that we have not achieved the required probability. It is, however, feasible that the default use of batching during training is interfering with estimating distributional properties of the data. Thus, we shall try again without batching.<jupyter_code>rbm = SequentialBernoulliRBM(
    num_output, num_input, 
    batch_size=1.0
)
rbm.fit(X_train)

preds = rbm.reconstruct(X_states)
probs = preds[:, 0]  # p(x_1)
probs = np.array(list(zip(1 - probs, probs)))
print(probs[range(2), X_states[:, 0]])<jupyter_output>[0.4721875 0.5278125]<jupyter_text>Here we have used the standard model which uses mean-field approximations to bypass the usual intractability of Boltzmann machine calculations (caused by the necessity of marginalising over many variables).As an alternative, let us now derive a sequence model that explicitly marginalises over the unknown output.We shall start from first principles, assuming the simplest possible model.  One-input, one-output To see what the exact gradient should be, let us start with the simplest Bernoulli RBM model, namely\begin{eqnarray}p(x_1,y_1) & = & \frac{e^{f(x_1,y_1)}}{\sum_{x_1'=0}^{1}\sum_{y_1'=0}^{1} e^{f(x_1',y_1')}}\,,\end{eqnarray}with\begin{eqnarray}f(x_1,y_1) & = & a_1 x_1 + x_1 W_{11} y_1 + b_1 y_1\,.\end{eqnarray}Thus\begin{eqnarray}p(x_1) & = & \frac{\sum_{y_1=0}^{1} e^{f(x_1,y_1)}}{\sum_{x_1'=0}^{1}\sum_{y_1'=0}^{1} e^{f(x_1',y_1')}}\\& = &\frac{e^{a_1 x_1}\left[1+e^{b_1 + x_1 W_{11}}\right]}{\sum_{x_1'=0}^{1}e^{a_1 x_1'}\left[1+e^{b_1 + x_1' W_{11}}\right]}\\& = &\frac{e^{a_1 x_1}\left[1+e^{b_1 + x_1 W_{11}}\right]}{\left[1+e^{b_1}\right]+e^{a_1}\left[1+e^{b_1 + W_{11}}\right]}\\& = &\frac{ e^{a_1 x_1+\ln\left[1+e^{b_1 + x_1 W_{11}}\right]-\ln\left[1+e^{b_1}\right]}}{ 1+e^{a_1+\ln\left[1+e^{b_1 + W_{11}}\right]-\ln\left[1+e^{b_1}\right]}}\,,\end{eqnarray}with the consequence that\begin{eqnarray}p(x_1=1) & = & \sigma\left(a_1+\ln\left[1+e^{b_1 + W_{11}}\right]-\ln\left[1+e^{b_1}\right]\right)\,,\end{eqnarray}where $\sigma(\cdot)$ is the logistic sigmoid function. Now let $p_1\doteq p(x_1=1)$, such that the log-likelihood of a single observation of $x_1$ is\begin{eqnarray}L & = & \ln\left[p_1^{x_1}\,(1-p_1)^{1-x_1}\right]\\& = &x_1\ln p_1 + (1-x_1)\ln(1-p_1)\,,\\\Rightarrow \nabla L & = &\frac{x_1}{p_1}\nabla p_1-\frac{1-x_1}{1-p_1}\nabla p_1\\& = &\frac{x-p_1}{p_1 (1-p_1)}\nabla p_1\,.\end{eqnarray} However, we recall that\begin{eqnarray}\nabla\sigma(\theta) & = & \sigma(\theta)\,[1-\sigma(\theta)]\,\nabla\theta\,,\\\Rightarrow\nabla p_1 & = & p_1 (1-p_1)\nabla\theta\,,\\\Rightarrow\nabla L & = & (x_1-p_1)\nabla\theta\,,\end{eqnarray}where\begin{eqnarray}\theta & = & a_1+\ln\left[1+e^{b_1 + W_{11}}\right]-\ln\left[1+e^{b_1}\right]\,.\end{eqnarray} Clearly, we have\begin{eqnarray}\frac{\partial\theta}{\partial a_1} & = & 1\,,\\\Rightarrow \frac{\partial L}{\partial a_1} & = & x_1 - p_1\,,\\\Rightarrow \frac{\partial\langle L\rangle}{\partial a_1} & = &\left\langle\frac{\partial L}{\partial a_1}\right\rangle = \langle x_1\rangle - p_1\,,\end{eqnarray}and thus the estimates of $a_1$ should converge when $p_1=\langle x_1\rangle$, which is just the empirical probability of the input bit being 1. Similarly, we have\begin{eqnarray}\frac{\partial\theta}{\partial W_{11}} & = &\frac{e^{b_1+W_{11}}}{1+e^{b_1+W_{11}}} = \sigma(b_1+W_{11})\,,\\\Rightarrow\frac{\partial L}{\partial W_{11}}& = & \left(x_1-p_1\right)\,\sigma(b_1+W_{11})\,,\end{eqnarray}and\begin{eqnarray}\frac{\partial\theta}{\partial b_{1}} & = &\frac{e^{b_1+W_{11}}}{1+e^{b_1+W_{11}}}-\frac{e^{b_1}}{1+e^{b_1}}\\& = &\sigma(b_1+W_{11})-\sigma(b_1)\,,\\\Rightarrow\frac{\partial L}{\partial b_{1}}& = & \left(x_1-p_1\right)\,\left[\sigma(b_1+W_{11})-\sigma(b_1)\right]\,.\end{eqnarray} Alternatively, we can go back to the model and examine the estimation approach to the gradients. Thus\begin{eqnarray}L & = & \ln\sum_{y_1=0}^{1} e^{f(x_1,y_1)}-\ln\sum_{x_1'=0}^{1}\sum_{y_1'=0}^{1} e^{f(x_1',y_1')}\,,\\\Rightarrow\nabla L & = & \mathbb{E}_{y_1\mid x_1}\left[\nabla f(x_1,y_1)\right]-\mathbb{E}_{x_1',y_1'}\left[\nabla f(x_1',y_1')\right]\,.\end{eqnarray} Now\begin{eqnarray}\frac{\partial f}{\partial a_1} & = & x_1\,,\\\Rightarrow\frac{\partial L}{\partial a_1} & = & \mathbb{E}_{y_1\mid x_1}\left[x_1\right]-\mathbb{E}_{x_1',y_1'}\left[x_1'\right]\\& = &x_1-p(x_1=1) = x_1-p_1\,.\end{eqnarray}In comparison, the mean-field gradient is\begin{eqnarray}\frac{\partial L}{\partial a_1} & = & x_1-\bar{x}_1'\,,\end{eqnarray}where\begin{eqnarray}\bar{x}_1' \doteq p(x_1=1\mid\bar{y}_1)\,, &&\bar{y}_1 \doteq p(y_1=1\mid x_1)\,.\end{eqnarray}Hence,we see that the reconstruction probability $\bar{x}_1'$ is now replaced by theexact probability $p_1$, denoted as $\bar{x}_1'\hookleftarrow p_1$. Similarly, we have\begin{eqnarray}\frac{\partial f}{\partial b_1} & = & y_1\,,\\\Rightarrow\frac{\partial L}{\partial b_1} & = & \mathbb{E}_{y_1\mid x_1}\left[y_1\right]-\mathbb{E}_{x_1',y_1'}\left[y_1'\right]\\& = & p(y_1=1\mid x_1)-p(y_1=1)\\& = & \bar{y}_1-q_1\,,\end{eqnarray}using the definition of $\bar{y}_1$ above,where we have now defined $q_1\doteq p(y_1=1)$ for reasons of symmetry.We deduce that\begin{eqnarray}p(y_1\mid x_1) & = & \frac{e^{f(x_1,y_1)}}{\sum_{y_1'=0}^{1} e^{f(x_1,y_1')}}\\& = &\frac{e^{a_1 x_1+x_1 W_{11}y_1+b_1 y_1}}{\sum_{y_1'=0}^{1} e^{a_1 x_1+x_1 W_{11}y_1'+b_1 y_1'}}\\& = &\frac{e^{(b_1+x_1 W_{11})y_1}}{1+e^{b_1+x_1 W_{11}}}\,,\\\Rightarrowp(y_1=1\mid x_1) & = & \sigma(b_1+x_1 W_{11})\,,\end{eqnarray}and therefore\begin{eqnarray}p(y_1=1) & = & p(y_1=1\mid x_1=1)\,p(x_1=1)+p(y_1=1\mid x_1=0)\,p(x_1=0)\\& = &p_1\sigma(b_1+W_{11})+(1-p_1)\sigma(b_1)\,,\end{eqnarray}giving\begin{eqnarray}\frac{\partial L}{\partial b_1} & = & \sigma(b_1+x_1 W_{11})-\left[p_1\sigma(b_1+W_{11})+(1-p_1)\sigma(b_1)\right]\\& = &(x_1-p_1)\left[\sigma(b_1+W_{11})-\sigma(b_1)\right]\,,\end{eqnarray}as expected. We note that the mean-field approximation is\begin{eqnarray}\frac{\partial L}{\partial b_1} & = & \bar{y}_1-\bar{y}_1'\,,\end{eqnarray}where\begin{eqnarray}\bar{y}_1' & \doteq & p(y_1=1\mid\bar{x}_1')\,.\end{eqnarray}Hence, the exact gradient is obtained via the replacement $\bar{y}_1'\hookleftarrow q_1$. Finally, we have\begin{eqnarray}\frac{\partial f}{\partial W_{11}} & = & x_1 y_1\,,\\\Rightarrow\frac{\partial L}{\partial W_{11}} & = & \mathbb{E}_{y_1\mid x_1}\left[x_1 y_1\right]-\mathbb{E}_{x_1',y_1'}\left[x_1' y_1'\right]\\& = &x_1\bar{y}_1-\mathbb{E}_{x_1'}\left[x_1' p(y_1'=1\mid x_1')\right]\\& = &x_1\bar{y}_1-p_1\,p(y_1=1\mid x_1=1)\,.\end{eqnarray}For completeness, we note that\begin{eqnarray}\frac{\partial L}{\partial W_{11}} & = & x_1\sigma(b_1+x_1 W_{11})-p_1\sigma(b_1+W_{11})\\& = & (x_1-p_1)\,\sigma(b_1+W_{11})\,,\end{eqnarray}as expected. In comparison, the mean-field approximation is \begin{eqnarray}\frac{\partial L}{\partial W_{11}} & = & x_1\bar{y}_1-\bar{x}_1'\bar{y}_1'\,, \end{eqnarray}with $\bar{y}_1$ and $\bar{y}_1'$ defined above.Thus, the exact gradient is obtained by the replacement$\bar{y}_1'\hookleftarrow p(y_1=1\mid x_1=1)=\sigma(b_1+W_{11})$. Note that thisis a different replacement for $\bar{y}_1'$ than that used for $\frac{\partial L}{\partial b_1}$!  Two-input, one-output Next, we extend the model such that the input is a sequence of two bits. We now have\begin{eqnarray}f(x_1,x_2,y_1) & = & a_1 x_1 + a_2 x_2 + x_1 W_{11} y_1 + x_2 W_{21} y_1 + b_1 y_1\,,\end{eqnarray}such that\begin{eqnarray}p(x_1,x_2,y_1) & = &\frac{e^{f(x_1,x_2,y_1)}}{\sum_{x_1'=0}^{1}\sum_{x_2'=0}^{1}\sum_{y_1'=0}^{1}e^{f(x_1',x_2',y_1')}}\,.\end{eqnarray}We note a property of the RBM that $x_1$ and $x_2$ are conditionally independent,namely\begin{eqnarray}p(x_1,x_2\mid y_1) & = &\frac{e^{f(x_1,x_2,y_1)}}{\sum_{x_1'=0}^{1}\sum_{x_2'=0}^{1}e^{f(x_1',x_2',y_1)}}\\& = & \frac{e^{(a_1 + W_{11}y_1) x_1+(a_2 + W_{21}y_1) x_2+b_1 y_1}}{\sum_{x_1'=0}^{1}\sum_{x_2'=0}^{1} e^{(a_1 + W_{11}y_1) x_1'+(a_2 + W_{21}y_1) x_2'+b_1 y_1}}\\& = &\frac{e^{(a_1 + W_{11}y_1) x_1}}{\sum_{x_1'=0}^{1}e^{(a_1 + W_{11}y_1) x_1'}}\frac{e^{(a_2 + W_{21}y_1) x_2}}{\sum_{x_2'=0}^{1}e^{(a_2 + W_{21}y_1) x_2'}}\\& = &p(x_1\mid y_1)\,p(x_2\mid y_1)\,,\end{eqnarray}such that\begin{eqnarray}p(x_i=1\mid y_1) & = & \sigma(a_i+W_{i1}y_1)\,.\end{eqnarray} However, $x_1$ and $x_2$ are **not** unconditionally independent, because\begin{eqnarray}p(x_2\mid x_1) & = & \frac{\sum_{y_1=0}^{1}e^{f(x_1,x_2,y_1)}}{\sum_{x_2'=0}^{1}\sum_{y_1'=0}^{1}e^{f(x_1,x_2',y_1')}}\\& = &\frac{\sum_{y_1=0}^{1}e^{a_1 x_1+a_2 x_2 + (b_1+x_1 W_{11}+x_2 W_{21}) y_1}}{\sum_{x_2'=0}^{1}\sum_{y_1'=0}^{1} e^{a_1 x_1+a_2 x_2' + (b_1+x_1 W_{11}+x_2' W_{21}) y_1'}}\\& = &\frac{e^{a_2 x_2}\sum_{y_1=0}^{1}e^{(b_1+x_1 W_{11}+x_2 W_{21}) y_1}}{\sum_{x_2'=0}^{1}e^{a_2 x_2'}\sum_{y_1'=0}^{1}e^{(b_1+x_1 W_{11}+x_2' W_{21}) y_1'}}\\& = &\frac{e^{a_2 x_2+\ln\left[1+e^{b_1+x_1 W_{11}+x_2 W_{21}}\right]}}{\sum_{x_2'=0}^{1}e^{a_2 x_2'+\ln\left[1+e^{b_1+x_1 W_{11}+x_2' W_{21}}\right]}}\\& = &\frac{e^{a_2 x_2+\ln\left[1+e^{b_1+x_1 W_{11}+x_2 W_{21}}\right]-\ln\left[1+e^{b_1+x_1 W_{11}}\right]}}{1+e^{a_2+\ln\left[1+e^{b_1+x_1 W_{11}+W_{21}}\right]-\ln\left[1+e^{b_1+x_1 W_{11}}\right]}}\,,\end{eqnarray}such that\begin{eqnarray}p(x_2=1\mid x_1) & = & \sigma\left(a_2+\ln\left[1+e^{b_1+x_1 W_{11}+W_{21}}\right]-\ln\left[1+e^{b_1+x_1 W_{11}}\right]\right)\,.\end{eqnarray} Now, due to this dependence, we note that the initial probability for the two-bit input model is in fact given by\begin{eqnarray}p(x_1=1) & = & \sum_{x_2=0}^{1} p(x_1=1,x_2)\,.\end{eqnarray}However, summing over unknown future values of a sequence is quickly going to become intractable for large sequences. Hence, we make an appeal to causality that future values cannot affect past values (although theywill affect *inferences* about past values). Thus, as an approximation, we assume the initial probabilityis just\begin{eqnarray}p(x_1=1) & \doteq & \sigma\left(a_1+\ln\left[1+e^{b_1+W_{11}}\right]-\ln\left[1+e^{b_1}\right]\right)\,,\end{eqnarray}as derived in the previous section for the one-bit input model.In effect, we truncate the full model of $f(\cdot)$ to only include as much of the input sequence as hasbeen observed.This causal approximation is the true basis of our sequential model.  $F$-input, $H$-output The full model, for input ${\bf x}\in\{0,1\}^{F}$ and output${\bf y}\in\{0,1\}^{H}$, is given by\begin{eqnarray}p({\bf x},{\bf y}) & = &\frac{e^{f({\bf x},{\bf y})}}{\sum_{{\bf x}'\in\{0,1\}^{F}}\sum_{{\bf y}'\in\{0,1\}^{H}}e^{f({\bf x}',{\bf y}')}}\,,\end{eqnarray}where\begin{eqnarray}f({\bf x},{\bf y}) & = & {\bf a}^{T}{\bf x}+{\bf x}^{T}{\bf W}{\bf y}+{\bf b}^{T}{\bf y}\,.\end{eqnarray} The sequential model for partial input ${\bf x}_{1:i}=(x_1,\ldots,x_i)$ istherefore (as a generalisation from the previous section) given by\begin{eqnarray}p({\bf x}_{1:i},{\bf y}) & = &\frac{e^{f({\bf x}_{1:i},{\bf y})}}{\sum_{{\bf x}_{1:i}'\in\{0,1\}^{i}}\sum_{{\bf y}'\in\{0,1\}^{H}}e^{f({\bf x}_{1:i}',{\bf y}')}}\,,\end{eqnarray}where we define the truncated model as\begin{eqnarray}f({\bf x}_{1:i},{\bf y}) & \doteq & {\bf a}_{1:i}^{T}{\bf x}_{1:i}+{\bf x}_{1:i}^{T}{\bf W}_{1:i,:}{\bf y}+{\bf b}^{T}{\bf y}\,.\end{eqnarray}Thus, we obtain\begin{eqnarray}p(x_i\mid{\bf x}_{1:i-1}) & = &\frac{\sum_{{\bf y}\in\{0,1\}^{H}}e^{f({\bf x}_{1:i-1},x_i,{\bf y})}}{\sum_{x_i'\in\{0,1\}}\sum_{{\bf y}'\in\{0,1\}^{H}}e^{f({\bf x}_{1:i-1},x_i',{\bf y}')}}\,,\\\Rightarrowp(x_i=1\mid{\bf x}_{1:i-1}) & = &\sigma\left( a_i+\sum_{j=1}^{H}\left\{  \ln\left[1+e^{B_{ij}+W_{ij}}\right]-\ln\left[1+e^{B_{ij}}\right] \right\}\right)\,,\end{eqnarray}where\begin{eqnarray}B_{ij} & \doteq & b_j+\sum_{k=1}^{i-1}x_k W_{kj}=b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}\,.\end{eqnarray} We now define the log-likelihood of the single sequence ${\bf x}$ as $L=\sum_{i=1}^{F} L_i$, for conditional log-likelihood\begin{eqnarray}L_i & \doteq & \ln p(x_i\mid{\bf x}_{1:i-1})\\& = &\ln\sum_{{\bf y}\in\{0,1\}^{H}}e^{f({\bf x}_{1:i-1},x_i,{\bf y})}-\ln\sum_{x_i'\in\{0,1\}}\sum_{{\bf y}'\in\{0,1\}^{H}}e^{f({\bf x}_{1:i-1},x_i',{\bf y}')}\,,\end{eqnarray}such that\begin{eqnarray}\nabla L & = & \sum_{i=1}^{F}\left\{\mathbb{E}_{{\bf y}\mid{\bf x}_{1:i}}\left[\nabla f({\bf x}_{1:i-1},x_i,{\bf y})\right]-\mathbb{E}_{x_i',{\bf y}'\mid{\bf x}_{1:i-1}}\left[  \nabla f({\bf x}_{1:i-1},x_i',{\bf y}')\right]\right\}\\& = &\sum_{i=1}^{F}\left\{\mathbb{E}_{{\bf y}\mid{\bf x}_{1:i-1},x_i}\left[  \nabla f({\bf x}_{1:i-1},x_i,{\bf y})\right]-\mathbb{E}_{x_i'\mid{\bf x}_{1:i-1}}\left[\mathbb{E}_{{\bf y}'\mid {\bf x}_{1:i-1},x_i'}\left[  \nabla f({\bf x}_{1:i-1},x_i',{\bf y}')\right]\right]\right\}\,.\end{eqnarray} We observe that for the $j$-th output variable $y_j$ we have\begin{eqnarray}\mathbb{E}_{{\bf y}\mid{\bf x}_{1:i}}[y_j] & = &p(y_j=1\mid{\bf x}_{1:i}) \doteq \bar{y}_j({\bf x}_{1:i-1},x_i)\,,\end{eqnarray}and\begin{eqnarray}\mathbb{E}_{x_i,{\bf y}\mid{\bf x}_{1:i-1}}\left[y_j\right]& = & p(y_j=1\mid{\bf x}_{1:i-1}) \doteq \tilde{y}_j({\bf x}_{1:i-1})\,,\end{eqnarray}where we shall derive the required formulae later. Thus, we obtain\begin{eqnarray}\frac{\partial f}{\partial b_j} = y_j & \Rightarrow &\frac{\partial L}{\partial b_j} = \sum_{i=1}^{F}\left\{\bar{y}_j({\bf x}_{1:i-1},x_i)-\tilde{y}_j({\bf x}_{1:i-1})\right\}\,.\end{eqnarray} However, due to our assumption of causality, we note that for the $k$-th input variable $x_k$ we have\begin{eqnarray}\mathbb{E}_{{\bf y}\mid{\bf x}_{1:i}}[x_k] & = &\left\{\begin{array}{ll}x_k\,, & k\le i\\0\,, & k>i\end{array}\right.\,,\end{eqnarray}and\begin{eqnarray}\mathbb{E}_{x_i\mid{\bf x}_{1:i-1}}[x_k] & = &\left\{\begin{array}{ll}x_k\,, & k<i\\p(x_i=1\mid{\bf x}_{1:i-1})\,, & k=i\\0\,, & k>i\end{array}\right.\,.\end{eqnarray}Thus, we obtain\begin{eqnarray}    \frac{\partial f}{\partial a_i} = x_i & \Rightarrow &\frac{\partial L}{\partial a_i} = x_i - \tilde{x}_i({\bf x}_{1:i-1})\,,\end{eqnarray}where we define\begin{eqnarray}    \tilde{x}_i({\bf x}_{1:i-1}) & \doteq &\mathbb{E}_{x_i\mid{\bf x}_{1:i-1}}[x_i] = p(x_i=1\mid{\bf x}_{1:i-1})\,.\end{eqnarray} We therefore deduce that\begin{eqnarray}\tilde{y}_j({\bf x}_{1:i-1}) & \doteq &\mathbb{E}_{x_i\mid{\bf x}_{1:i-1}}\left[\mathbb{E}_{{\bf y}\mid{\bf x}_{1:i-1},x_i}[y_j] \right]\\& = &\mathbb{E}_{x_i\mid{\bf x}_{1:i-1}}\left[p(y_j=1\mid{\bf x}_{1:i-1},x_i)\right]\\& = &p(x_i=0\mid{\bf x}_{1:i-1})\,p(y_j=1\mid{\bf x}_{1:i-1},x_i=0)+p(x_i=1\mid{\bf x}_{1:i-1})\,p(y_j=1\mid{\bf x}_{1:i-1},x_i=1)\\& = &[1-\tilde{x}_i({\bf x}_{1:i-1})]\,\bar{y}_j({\bf x}_{1:i-1},0)+\tilde{x}_i({\bf x}_{1:i-1})\,\bar{y}_j({\bf x}_{1:i-1},1)\,.\end{eqnarray} Lastly, we examine the interaction between the $k$-th input and the $j$-th output. We observe that\begin{eqnarray}\mathbb{E}_{{\bf y}\mid{\bf x}_{1:i}}[x_k y_j] & = &\left\{\begin{array}{ll}x_k\,\bar{y}_j({\bf x}_{1:i-1},x_i)\,, & k\le i\\0\,, & k>i\end{array}\right.\,,\end{eqnarray}and therefore\begin{eqnarray}\mathbb{E}_{x_i\mid{\bf x}_{1:i-1}}\left[\mathbb{E}_{{\bf y}\mid{\bf x}_{1:i}}[x_k y_j]\right] & = &\left\{\begin{array}{ll}x_k\,\tilde{y}_j({\bf x}_{1:i-1})\,, & k<i\\\tilde{x}_i({\bf x}_{1:i-1})\,\bar{y}_j({\bf x}_{1:i-1},1)\,, & k=i\\0\,, & k>i\end{array}\right.\,.\end{eqnarray} Consequently, we observe that\begin{eqnarray}\frac{\partial f}{\partial W_{kj}} & = &x_k\,y_j\\\Rightarrow\frac{\partial L}{\partial W_{kj}} & = &x_k\,\bar{y}_j({\bf x}_{1:k-1},x_k)-\tilde{x}_k({\bf x}_{1:k-1})\,\bar{y}_j({\bf x}_{1:k-1},1)+ x_k\sum_{i=k+1}^{F}\left[  \bar{y}_j({\bf x}_{1:i-1},x_i)-\tilde{y}_j({\bf x}_{1:i-1})\right]\,.\end{eqnarray}Swapping indices $i$ and $k$ for convenience, this gives the more familiar form of\begin{eqnarray}\frac{\partial L}{\partial W_{ij}} & = &x_i\,\bar{y}_j({\bf x}_{1:i-1},x_i)-\tilde{x}_i({\bf x}_{1:i-1})\,\bar{y}_j({\bf x}_{1:i-1},1)+ x_i\sum_{k=i+1}^{F}\left[  \bar{y}_j({\bf x}_{1:k-1},x_k)-\tilde{y}_j({\bf x}_{1:k-1})\right]\,.\end{eqnarray}Thus, we see that inference of $a_i$ relies only on past observations ${\bf x}_{1:i-1}$ and present observation $x_i$, but that inferences of $b_j$ and $W_{ij}$ also rely on future observations. Thisnicely demonstrates what we said earlier about causality, in that future observations cannot physically affectpast observations, but they can affect inferences about past observations. Finally, for completeness, we derive the missing formula for $\bar{y}_j({\bf x}_{1:i-1},x_i)$ by noting that\begin{eqnarray}p({\bf y}\mid{\bf x}_{1:i}) & = &\frac{e^{f({\bf x}_{1:i},{\bf y})}}{\sum_{{\bf y}'\in\{0,1\}^{H}}e^{f({\bf x}_{1:i},{\bf y}')}}\\& = &\frac{e^{ {\bf a}_{1:i}^{T}{\bf x}_{1:i}+{\bf x}_{1:i}^{T}{\bf W}_{1:i,:}{\bf y}+{\bf b}^{T}{\bf y}}}{\sum_{{\bf y}'\in\{0,1\}^{H}}e^{ {\bf a}_{1:i}^{T}{\bf x}_{1:i}+{\bf x}_{1:i}^{T}{\bf W}_{1:i,:}{\bf y}'+{\bf b}^{T}{\bf y}'}}\\& = &\frac{\prod_{j=1}^{H}e^{(b_j+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j})y_j}}{\prod_{j=1}^{H}\sum_{y_j=0}^{1}e^{(b_j+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j})y_j}}\\& = &\prod_{j=1}^{H}\frac{e^{(b_j+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j})y_j}}{1+e^{b_j+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j}}}\\\Rightarrowp(y_j=1\mid{\bf x}_{1:i})& = & \sigma\left(b_j+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j}\right)\\\Rightarrow\bar{y}_j({\bf x}_{1:i-1},x_i) & = & \sigma\left(B_{ij}+x_i W_{ij}\right)\,,\end{eqnarray}where $B_{ij}$ was defined earlier.  $F$-input, $H$-output Classifier We now consider the special case where the binary output vector ${\bf y}$ is restricted to being a$1$-of-$H$ or one-hot vector, such that ${\bf y}\in{\cal Y}_H$ where ${\cal Y}_H=\left\{{\bf y}\in\{0,1\}^{H}\mid\sum_{j=1}^{H}y_j=1\right\}$.This model corresponds to a $H$-class classifier. The predictive model is therefore\begin{eqnarray}p(x_i\mid{\bf x}_{1:i-1}) & = &\frac{\sum_{{\bf y}\in{\cal Y}_{H}}e^{f({\bf x}_{1:i-1},x_i,{\bf y})}}{\sum_{x_i'\in\{0,1\}}\sum_{{\bf y}'\in{\cal Y}_{H}}e^{f({\bf x}_{1:i-1},x_i',{\bf y}')}}\\& = &\frac{\sum_{{\bf y}\in{\cal Y}_{H}}e^{ {\bf a}_{1:i-1}^{T}{\bf x}_{1:i-1}+a_i x_i+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,:}{\bf y} +x_i W_{i,:}{\bf y}+{\bf b}^{T}{\bf y}}}{\sum_{x_i'\in\{0,1\}}\sum_{{\bf y}'\in{\cal Y}_{H}}e^{ {\bf a}_{1:i-1}^{T}{\bf x}_{1:i-1}+a_i x_i'+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,:}{\bf y}' +x_i' W_{i,:}{\bf y}'+{\bf b}^{T}{\bf y}'}}\\& = &\frac{ e^{a_i x_i} \sum_{{\bf y}\in{\cal Y}_{H}}e^{ \left({\bf b}^{T}+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,:}+x_i W_{i,:}\right){\bf y}}}{\sum_{x_i'\in\{0,1\}}e^{a_i x_i'}\sum_{{\bf y}'\in{\cal Y}_{H}}e^{ \left({\bf b}^{T}+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,:}+x_i' W_{i,:}\right){\bf y}'}}\\& = &\frac{ e^{a_i x_i} \sum_{j=1}^{H}e^{  b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}+x_i W_{ij}}}{\sum_{x_i'\in\{0,1\}}e^{a_i x_i'}\sum_{j=1}^{H}e^{  b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}+x_i' W_{ij}}}\\& = &\frac{ e^{a_i x_i+\ln \sum_{j=1}^{H}e^{  b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}+x_i W_{ij}}}}{ e^{\ln  \sum_{j=1}^{H}e^{   b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}  } }+ e^{a_i+\ln  \sum_{j=1}^{H}e^{   b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}+W_{ij}  } }}\\& = &\frac{ e^{a_i x_i  +\ln\sum_{j=1}^{H}e^{b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}+x_i W_{ij}}  -\ln\sum_{j=1}^{H}e^{b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}} }}{ 1+e^{  a_i  +\ln\sum_{j=1}^{H}e^{b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}+W_{ij}}  -\ln\sum_{j=1}^{H}e^{b_j+{\bf x}_{1:i-1}^{T}{\bf W}_{1:i-1,j}} }}\,.\end{eqnarray} Consequently, we have\begin{eqnarray}\tilde{x}_i({\bf x}_{1:i-1}) & \doteq & p(x_i=1\mid{\bf x}_{1:i-1})=\sigma\left(a_i+\ln\sum_{j=1}^{H}e^{B_{ij}+W_{i,j}}-\ln\sum_{j=1}^{H}e^{B_{ij}}\right)\,,\end{eqnarray}where $B_{ij}$ was defined in the previous section. Similarly, we have\begin{eqnarray}p({\bf y}\mid{\bf x}_{1:i}) & = &\frac{e^{f({\bf x}_{1:i},{\bf y})}}{\sum_{{\bf y}'\in{\cal Y}_{H}}e^{f({\bf x}_{1:i},{\bf y}')}}\\& = &\frac{e^{ {\bf a}_{1:i}^{T}{\bf x}_{1:i}+{\bf x}_{1:i}^{T}{\bf W}_{1:i,:}{\bf y}+{\bf b}^{T}{\bf y}}}{\sum_{{\bf y}'\in{\cal Y}_{H}}e^{ {\bf a}_{1:i}^{T}{\bf x}_{1:i}+{\bf x}_{1:i}^{T}{\bf W}_{1:i,:}{\bf y}'+{\bf b}^{T}{\bf y}'}}\\& = &\frac{\prod_{j=1}^{H}e^{(b_j+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j})y_j}}{\sum_{j=1}^{H}e^{b_j+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j}}}\\\Rightarrowp(y_j=1\mid{\bf x}_{1:i})& = & \frac{e^{b_j+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j}}}{\sum_{j'=1}^{H}e^{b_{j'}+{\bf x}_{1:i}^{T}{\bf W}_{1:i,j'}}}\\\Rightarrow\bar{y}_j({\bf x}_{1:i-1},x_i) & = & \frac{e^{B_{ij}+x_i W_{ij}}}{\sum_{j'=1}^{H}e^{B_{ij'}+x_i W_{ij'}}}\,.\end{eqnarray}  Continue the Experiment Now let us continue the experiment by training on the exact model (which sums over all possible outputs ${\bf y}$), rather than using the mean-field approximation.<jupyter_code>rbm = ExactSequentialBernoulliRBM(
    num_output, num_input, 
    batch_size=1.0
)
rbm.fit(X_train)

preds = rbm.reconstruct(X_states)
probs = preds[:, 0]  # p(x_1)
probs = np.array(list(zip(1 - probs, probs)))
print(probs[range(2), X_states[:, 0]])
rbm._compute_score_gradients(X_train)<jupyter_output><empty_output><jupyter_text>We see that the model appears to have converged, and we obtained the empirical probability.<jupyter_code>print("p(x_1=1) = %f" % p1)<jupyter_output>p(x_1=1) = 0.527813<jupyter_text>2-input, 1-output Let us now examine the difference, if any, between the two models (exact and mean-field) for the case of 2 binary inputs. We retain 1 binary output for simplicity.<jupyter_code>num_output = 1
num_input = 2
X_train = _window(b_series, num_input)
X_states = np.array([[0,0],[0,1], [1,0], [1,1]])
ind1 = X_train[:, 0] == 1
p1 = sum(ind1) / len(ind1)
ind2 = X_train[:, 1] == 1
p2 = sum(ind2) / len(ind2)
print("p(x_1=1) = %f, p(x_2=1) = %f" % (p1, p2))
p2g0 = sum(ind2 & ~ind1) / sum(~ind1) 
p2g1 = sum(ind2 & ind1) / sum(ind1)
print("p(x_2=1|x_1=0) = %f, p(x_2=1|x_1=1) = %f" % (p2g0, p2g1))
rbm = ExactSequentialBernoulliRBM(
    num_output, num_input,
    n_iter=10000,
    batch_size=1.0
)
rbm.fit(X_train)

preds = rbm.reconstruct(X_states)
print(preds)
rbm._compute_score_gradients(X_train)<jupyter_output><empty_output><jupyter_text>Note that we require about 10,000 iterations to converge to the empirical values, even for this simple model.Clearly, some form of gradient acceleration would be useful.<jupyter_code>rbm = SequentialBernoulliRBM(
    num_output, num_input,
    n_iter=10000,
    batch_size=1.0
)
rbm.fit(X_train)

preds = rbm.reconstruct(X_states)
print(preds)
rbm._compute_score_gradients(X_train)<jupyter_output><empty_output><jupyter_text>Surprisingly, there doesn't seem to be any difference between the two models. And now for something completely different - let's take a quick look at gradient acceleration.In particular, we use a 1-step LBFGS algorithm with optional line-search (to ensure the `x_score` always increases). To explain the algorithm, we use [Wikipedia](https://en.wikipedia.org/wiki/Limited-memory_BFGS) as a reference, includingcopying its notation, but assume exactly $m=1$ step.Imagine we have some parameter estimate ${\bf x}_{k-1}$, andthen move to a new estimate ${\bf x}_k$. The second-order Taylor series expansion is then\begin{eqnarray}f({\bf x}_k) & = & f({\bf x}_{k-1})+\nabla^{T}f({\bf x}_{k-1})\,({\bf x}_k-{\bf x}_{k-1})+({\bf x}_k-{\bf x}_{k-1})^{T}\nabla\nabla^{T}f({\bf x}_{k-1})\,({\bf x}_k-{\bf x}_{k-1})+O(\|{\bf x}_k-{\bf x}_{k-1}\|^3)\\& = & f({\bf x}_{k-1})+{\bf g}_{k-1}^{T}{\bf s}_{k-1}+{\bf s}_{k-1}^{T}{\bf H}_{k-1}^{-1}{\bf s}_{k-1}+O(\|{\bf s}_{k-1}\|^3)\,,\end{eqnarray}and its gradient is\begin{eqnarray}\nabla f({\bf x}_k) & = & \nabla f({\bf x}_{k-1})+\nabla\nabla^{T}f({\bf x}_{k-1})\,({\bf x}_k-{\bf x}_{k-1})+O(\|{\bf x}_k-{\bf x}_{k-1}\|^2)\\\Rightarrow{\bf g}_k & = &{\bf g}_{k-1}+{\bf H}_{k-1}^{-1}{\bf s}_{k-1}+O(\|{\bf s}_{k-1}\|^2)\,.\end{eqnarray}Letting ${\bf y}_{k-1}\doteq{\bf g}_k-{\bf g}_{k-1}$ then gives\begin{eqnarray}{\bf y}_{k-1} \approx {\bf H}_{k-1}^{-1}{\bf s}_{k-1}& \Rightarrow & {\bf y}_{k-1}^{T}{\bf H}_{k-1}{\bf y}_{k-1} \approx {\bf y}_{k-1}^{T}{\bf s}_{k-1}\,.\end{eqnarray}For convenience, we assume\begin{eqnarray}{\bf H}_{k-1}\approx\gamma_k{\bf I} & \Rightarrow &\gamma_k \doteq \frac{{\bf y}_{k-1}^{T}{\bf s}_{k-1}}{{\bf y}_{k-1}^{T}{\bf y}_{k-1}}\,.\end{eqnarray}Note that [Wikipedia](https://en.wikipedia.org/wiki/Limited-memory_BFGS) called this${\bf H}^0_k$ rather than ${\bf H}_{k-1}$, in order to allow for further (i.e. $m>1$)steps into the past; presumably this is why they defined $\gamma_k$ rather than $\gamma_{k-1}$. We now seek an update to the next parameter estimate ${\bf x}_{k+1}$.This is obtained from the gradient equation\begin{eqnarray}\nabla f({\bf x}_{k+1}) & = & \nabla f({\bf x}_{k})+\nabla\nabla^{T}f({\bf x}_{k})\,({\bf x}_{k+1}-{\bf x}_{k})+O(\|{\bf x}_{k+1}-{\bf x}_{k}\|^2)\end{eqnarray}by assuming that $\nabla f({\bf x}_{k+1})\approx{\bf 0}$,giving\begin{eqnarray}{\bf x}_{k+1} & \approx & {\bf x}_k-\left[\nabla\nabla^T f({\bf x}_k)\right]^{-1}\,\nabla f({\bf x}_k)={\bf x}_k-{\bf H}_k{\bf g}_k\,.\end{eqnarray}Thus, the Newton-Raphson update direction is just ${\bf d}_k=-{\bf H}_k{\bf g}_k$. Note that the actual update will typically be some${\bf s}_k=\epsilon\,{\bf d}_k$, with step-size $\epsilon$, to allow for line searching.  The BFGS update from ${\bf H}_{k-1}$ to ${\bf H}_k$ is given by\begin{eqnarray}{\bf H}_k & = &\left(I-\rho_{k-1}{\bf s}_{k-1}{\bf y}_{k-1}^T\right){\bf H}_{k-1}\left(I-\rho_{k-1}{\bf y}_{k-1}{\bf s}_{k-1}^T\right)+\rho_{k-1}{\bf s}_{k-1}{\bf s}_{k-1}^{T}\,,\end{eqnarray}where $\rho_{k-1}\doteq\frac{1}{{\bf y}_{k-1}^T{\bf s}_{k-1}}$.Thus, the update direction is\begin{eqnarray}{\bf d}_k & = & -{\bf H}_k{\bf g}_k\\& \approx &-\left(I-\rho_{k-1}{\bf s}_{k-1}{\bf y}_{k-1}^T\right)\gamma_{k}{\bf I}\left(I-\rho_{k-1}{\bf y}_{k-1}{\bf s}_{k-1}^T\right){\bf g}_k-\rho_{k-1}{\bf s}_{k-1}{\bf s}_{k-1}^{T}{\bf g}_k\,.\end{eqnarray}For convenience, we define $\alpha_{k-1}\doteq\rho_{k-1}{\bf s}_{k-1}^{T}{\bf g}_k$,and\begin{eqnarray}{\bf q}_{k-1} & \doteq & \left(I-\rho_{k-1}{\bf y}_{k-1}{\bf s}_{k-1}^T\right){\bf g}_k={\bf g}_k-\alpha_{k-1}{\bf y}_{k-1}\,,\end{eqnarray}resulting in\begin{eqnarray}{\bf d}_k & \approx &-\left(I-\rho_{k-1}{\bf s}_{k-1}{\bf y}_{k-1}^T\right)\gamma_{k}{\bf q}_{k-1}-\alpha_{k-1}{\bf s}_{k-1}\,.\end{eqnarray}Similarly, we define $\beta_{k-1}\doteq\gamma_k\rho_{k-1}{\bf y}_{k-1}^T{\bf q}_{k-1}$ for convenience,giving\begin{eqnarray}{\bf d}_k & \approx &-\gamma_{k}{\bf q}_{k-1}+\beta_{k-1}{\bf s}_{k-1}-\alpha_{k-1}{\bf s}_{k-1}\,.\end{eqnarray}Note that this direction is called ${\bf z}$ in [Wikipedia](https://en.wikipedia.org/wiki/Limited-memory_BFGS).<jupyter_code>def diff(a, b):
    return tuple(_a - _b for _a, _b in zip(a, b))
def dot(a, b):
    return sum(np.vdot(_a, _b) for _a, _b in zip(a, b))
def mult(s, v):
    return tuple(s * _v for _v in v)
def add(a, b):
    return tuple(_a + _b for _a, _b in zip(a, b))
def grad_accel(s0, y0, g1):
    #s0 = diff(x1, x0)
    #y0 = diff(g1, g0)
    k1 = dot(s0, y0)  # 1 / rho_0
    k2 = dot(y0, y0)
    k3 = dot(s0, g1)
    alpha_0 = k3 / k1
    gamma_1 = k1 / k2
    d = mult(-gamma_1, diff(g1, mult(alpha_0, y0)))
    beta0 = -dot(y0, d) / k1
    d = add(d, mult(beta0 - alpha_0, s0))
    return d
def max_abs(x):
    return np.max([np.max(np.abs(_x)) for _x in x])
rbm = SequentialBernoulliRBM(
    num_output, num_input,
    n_iter=1,
    batch_size=1.0
)
num_iters = 0
score0 = rbm.score(X_train)['x_score']
print("iter=%2d, score=%10.8f" % (num_iters, score0))
x0 = rbm.copy_parameters()
g0 = rbm._compute_score_gradients(X_train)
rho0 = rbm.step_size
num_iters += 1
# Perform simple gradient ascent
x1 = add(x0, mult(rho0, g0))
num_steps = 1
rbm._set_parameters(*x1)
score1 = rbm.score(X_train)['x_score']
print("iter=%2d, score=%10.8f" % (num_iters, score1))
g1 = rbm._compute_score_gradients(X_train)
while True:
    s0 = diff(x1, x0)
    y0 = diff(g1, g0)
    if max_abs(y0) < 1e-8: # Any smaller leads to instability
        # Gradient has barely changed
        #   - stop to prevent division by zero
        break
    # Compute quasi-Newton gradient direction
    d = grad_accel(s0, y0, g1)
    x0 = x1
    g0 = g1
    score0 = score1
    ls = False
    while True:
        # Perform optional line-search until objective score improves
        x1 = add(x0, d)
        num_steps += 1
        rbm._set_parameters(*x1)
        score1 = rbm.score(X_train)['x_score']
        if score1 < score0:
            ls = True
            d = mult(0.5, d)
        else:
            break
    g1 = rbm._compute_score_gradients(X_train)
    num_iters += 1
    print("iter=%2d, score=%10.8f%s" % (num_iters, score1, "*" if ls else ""))
print("total steps=%d" % num_steps)
rbm._compute_score_gradients(X_train)
preds = rbm.reconstruct(X_states)
print(preds)<jupyter_output>[[0.52797749 0.55695367]
 [0.52797749 0.55695367]
 [0.52797749 0.50148015]
 [0.52797749 0.50148015]]<jupyter_text>Although the number of iterations required depends strongly on the random initialisation,we see that even 1-step LBFGS is significantly better than the (approximately) 10,000 iterations required without gradient acceleration.  We note that gradient acceleration schemes tend to exacerbate overshooting, which is manifested by the log-likelihood score oscillating rather than smoothly increasing.To prevent this, we require knowledge of the score being optimised, in order to performline search along the quasi-gradient direction. Ordinary gradient ascent can also overshoot, but typically suffers from undershooting.For our RBMs, this means knowing in advance whether we are optimising `x_score`, `y_score` or `xy_score`, which requires more detailed information about the RBM settings and computations.Alternatively, we could potentially detect oscillation without knowledge of the score, in the case where consecutive gradients point in "opposite" directions. Let's try a simple acceleration scheme.<jupyter_code>rbm = SequentialBernoulliRBM(
    num_output, num_input,
    n_iter=1,
    batch_size=1.0
)

rho0 = rbm.step_size

num_iters = 0
num_steps = 0
num_grads = 0

score = rbm.score(X_train)['x_score']  # not used in algorithm
print("DEBUG: iter=%3d, score=%10.8f, rho=%f" % (num_iters, score, rho0))

x0 = rbm.copy_parameters()
g0 = rbm._compute_score_gradients(X_train)
num_grads += 1
g0_g0 = dot(g0, g0)
while np.sqrt(g0_g0) > 1e-8:
    rho = rho0
    ls1 = False
    while True:
        # Simple gradient step
        x1 = add(x0, mult(rho, g0))
        rbm._set_parameters(*x1)
        num_steps += 1
        g1 = rbm._compute_score_gradients(X_train)
        num_grads += 1
        # Estimate quadratic turning point
        g1_g0 = dot(g1, g0)
        rho_star = -rho0 * g0_g0 / (g1_g0 - g0_g0)
        if rho_star < 0:
            # Overshoot - perform line search
            rho *= 0.5
            ls1 = True
        else:
            rho = rho_star
            break
    ls2 = False
    while True:
        # Accelerated gradient step
        x1 = add(x0, mult(rho, g0))
        rbm._set_parameters(*x1)
        num_steps += 1
        g1 = rbm._compute_score_gradients(X_train)
        num_grads += 1
        g1_g0 = dot(g1, g0)
        if g1_g0 < 0:
            # Oscillation - perform line search
            rho *= 0.5
            ls2 = True
        else:
            break
    # Reset at new position
    num_iters += 1
    ls1 = ("-" if ls1 else "")
    ls2 = ("+" if ls2 else "")
    score = rbm.score(X_train)['x_score']  # not used in algorithm
    print(
        "DEBUG: iter=%3d, score=%10.8f, rho=%f%s%s" 
        % (num_iters, score, rho, ls1, ls2)
    )
    x0 = x1
    g0 = g1
    g0_g0 = dot(g0, g0)
print("num_iters=%d, num_steps=%d, num_grads=%d" % (num_iters, num_steps, num_grads))
rbm._compute_score_gradients(X_train)
preds = rbm.reconstruct(X_states)
print(preds)<jupyter_output>[[0.5279775  0.55695357]
 [0.5279775  0.55695357]
 [0.5279775  0.50148022]
 [0.5279775  0.50148022]]<jupyter_text>We see that even this simplistic acceleration, which computes at least two gradients per iteration, greatly improves on having no acceleration.Note that we also seem to have prevented score oscillation, even without the algorithm actually computing the score. In fact, computing the score is a distraction, since it appears to converge way before the gradient shrinks enough to halt. Remember that it is the gradient that controls how accurately the model predictions match the empirical probabilities. Of course, this also essentially means we are deliberately (over?) fitting the training data, so the actual convergence test must remain domain-sepcific.  Future prediction Let us now consider a practical prediction task. We first train the RBM model on fortnightly binary sequences of daily stock rises and/or falls. We then test the model by using weekly sequences to predict the characteristics of the following weeks. For convenience, we compute the proportion of daily rises in the second week as a coarse measure of the weekly trend.<jupyter_code>num_output = 7
num_input = 14
X_train = _window(b_series, num_input)
rbm = SequentialBernoulliRBM(
    num_output, num_input,
    n_iter=1,
    batch_size=1.0
)
num_iters = 0
score0 = rbm.score(X_train)['x_score']
print("iter=%2d, score=%10.8f" % (num_iters, score0))
x0 = rbm.copy_parameters()
g0 = rbm._compute_score_gradients(X_train)
rho0 = rbm.step_size
num_iters += 1
# Perform simple gradient ascent
x1 = add(x0, mult(rho0, g0))
num_steps = 1
rbm._set_parameters(*x1)
score1 = rbm.score(X_train)['x_score']
print("iter=%2d, score=%10.8f" % (num_iters, score1))
g1 = rbm._compute_score_gradients(X_train)
while True:
    s0 = diff(x1, x0)
    y0 = diff(g1, g0)
    if max_abs(y0) < 1e-8: # Any smaller leads to instability
        # Gradient has barely changed
        #   - stop to prevent division by zero
        break
    # Compute quasi-Newton gradient direction
    d = grad_accel(s0, y0, g1)
    x0 = x1
    g0 = g1
    score0 = score1
    ls = False
    while True:
        # Perform optional line-search until objective score improves
        x1 = add(x0, d)
        num_steps += 1
        rbm._set_parameters(*x1)
        score1 = rbm.score(X_train)['x_score']
        if score1 < score0:
            ls = True
            d = mult(0.5, d)
        else:
            break
    g1 = rbm._compute_score_gradients(X_train)
    num_iters += 1
    if num_iters % 1000 == 0:
        print("iter=%2d, score=%10.8f%s" % (num_iters, score1, "*" if ls else ""))
print("total steps=%d" % num_steps)
P_train = np.mean(X_train[:, num_input // 2 :], axis=1)
X_test = X_train[:, 0 : num_input // 2]
X_pred = rbm.reconstruct(X_test)
P_test = np.mean(X_pred[:, num_input // 2 :], axis=1)
X_pred[0:10,:]
P_train[0:10]
X_train[0:10,:]
P_test[0:10]
acc = ((P_train <= 0.5) & (P_test <= 0.5)) | ((P_train > 0.5) & (P_test > 0.5))
np.mean(acc)<jupyter_output><empty_output><jupyter_text>We now compare this to the empirical baseline.<jupyter_code>P_train0 = np.mean(X_train[:, : num_input // 2], axis=1)
count00 = np.sum((P_train0 <= 0.5) & (P_train <= 0.5))
count01 = np.sum((P_train0 <= 0.5) & (P_train > 0.5))
count10 = np.sum((P_train0 > 0.5) & (P_train <= 0.5))
count11 = np.sum((P_train0 > 0.5) & (P_train > 0.5))
A = np.array([[count00, count01], [count10, count11]])
print(A)
print(np.sum(A, axis=1))
P = np.array([
    count01 / (count00 + count01),
    count11 / (count10 + count11)
])
print(P)
np.mean(P_train0 <= 0.5) * P[0] + np.mean(P_train0 > 0.5) * P[1]
np.sum(A, axis=0) / np.sum(A)<jupyter_output><empty_output><jupyter_text>Unfortunately, the predictive model appears to be less accurate than the baseline!  Sequence clustering Another possible use for the sequential model is clustering. What patterns can we discover amongst the various fortnighly sequences? Firstly, we shall examine the various scores.<jupyter_code>rbm = SequentialBernoulliRBM(
    num_output, num_input,
    n_iter=500, n_report=100,
    batch_size=1.0,
    is_RBC=True
)
rbm.fit(X_train)<jupyter_output>Iteration 0: xy_score=-11.650081, y_score=-1.945817, x_score=-9.704264, rmse=1.870856, mae=7.168183
Iteration 100: xy_score=-11.627625, y_score=-1.945837, x_score=-9.681788, rmse=1.867852, mae=6.605271
Iteration 200: xy_score=-11.627635, y_score=-1.945849, x_score=-9.681786, rmse=1.867851, mae=6.605271
Iteration 300: xy_score=-11.627630, y_score=-1.945846, x_score=-9.681784, rmse=1.867851, mae=6.605271
Iteration 400: xy_score=-11.627607, y_score=-1.945827, x_score=-9.681780, rmse=1.867851, mae=6.605271
Iteration 500: xy_score=-11.627562, y_score=-1.945788, x_score=-9.681774, rmse=1.867850, mae=6.605271<jupyter_text>Although the implementation maximises $p({\bf x})$, it appears that $p({\bf y}\mid{\bf x})$ might also be maximised (after a bit of struggling).<jupyter_code>rbm.n_iter=10000
rbm.n_report=1000
rbm.fit(X_train)<jupyter_output>Iteration 0: xy_score=-11.627562, y_score=-1.945788, x_score=-9.681774, rmse=1.867850, mae=6.605271
Iteration 1000: xy_score=-11.595553, y_score=-1.917486, x_score=-9.678066, rmse=1.867355, mae=6.605271
Iteration 2000: xy_score=-11.334338, y_score=-1.671932, x_score=-9.662407, rmse=1.865266, mae=6.567305
Iteration 3000: xy_score=-11.221774, y_score=-1.565064, x_score=-9.656711, rmse=1.864511, mae=6.555067
Iteration 4000: xy_score=-11.172699, y_score=-1.519828, x_score=-9.652871, rmse=1.864005, mae=6.536241
Iteration 5000: xy_score=-11.131369, y_score=-1.484020, x_score=-9.647348, rmse=1.863275, mae=6.532789
Iteration 6000: xy_score=-11.066540, y_score=-1.423568, x_score=-9.642971, rmse=1.862700, mae=6.514277
Iteration 7000: xy_score=-11.038727, y_score=-1.396842, x_score=-9.641886, rmse=1.862560, mae=6.516159
Iteration 8000: xy_score=-11.025159, y_score=-1.383549, x_score=-9.641611, rmse=1.862526, mae=6.534672
Iteration 9000: xy_score=-11.013104, y_score=-1.371242, x_score=-9.641862, r[...]<jupyter_text>Hmm, the `x_score` has started decreasing again! Also note that other experiments (not shown here) seem to indicate that gradient acceleration doesn't work well for the sequence classifier.<jupyter_code>Y_preds = rbm.predict(X_train)
Y_clust = np.argmax(Y_preds, axis=1)
Y_preds[0:10]
Y_clust[0:10]
X_train[0:10,:]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.23175702219603223, 'fraction_numerical': 0.0606462384600275, 'mean_word_length': 4.606965859030837, 'pattern_counts': {'":': 0, '<': 72, '<?xml version=': 0, '>': 71, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '420887', 'license': 'MIT', 'n_tokens_mistral': 18098, 'n_tokens_neox': 17160, 'n_words': 3378, 'path': 'RBMModels/notebooks/sequential_rbm.ipynb', 'repo_name': 'gaj67/gaj-data-science', 'size': '40728'}