<jupyter_start><jupyter_text>LicenseCopyright 2020 <NAME>Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License athttp://www.apache.org/licenses/LICENSE-2.0Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.  Automatic feature engineering using deep learning and Bayesian inference Author: <NAME>We will explore the use of autoencoders for automatic feature engineering. The idea is to automatically learn a set of features from raw data in an unsupervised manner that can be useful later in supervised learning tasks. Unlabelled data is inexpensive to obtain in large quantities therefore training a feature extractor in an unsupervised manner is a powerful first step towards later training supervised learning models which perhaps may not have access to a large amount of labelled data. Computer VisionWe will use the MNIST dataset for this purpose where the raw data is a 2 dimensional tensor of pixel intensities per image. The image is our unit of analysis: We will predict the probability of each class for each image. This is a multiclass classification task and we will use the accuracy score to assess model performance on the test fold.![](pixel_lattice.png)<jupyter_code>import inspect

import numpy as np

np.set_printoptions(suppress=True)

import pandas as pd

import torch
import torchvision

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")

from sklearn import linear_model
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import MinMaxScaler
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE

from scipy.stats import norm

import matplotlib.pyplot as plt
import plotnine

import typeguard

print(f"PyTorch version: {torch.__version__}")
print(f"Device: {device}")

data_path = "/home/hamaad/Projects/autoencoders_pytorch/data"
trainset = torchvision.datasets.MNIST(
    root=data_path,
    train=True,
    download=True,
    transform=torchvision.transforms.Compose(
        transforms=[torchvision.transforms.ToTensor()]
    ),
)

testset = torchvision.datasets.MNIST(
    root=data_path,
    train=False,
    download=True,
    transform=torchvision.transforms.Compose(
        transforms=[torchvision.transforms.ToTensor()]
    ),
)

batch_size = 100
trainloader = torch.utils.data.DataLoader(
    dataset=trainset, batch_size=batch_size, shuffle=True, num_workers=1
)
testloader = torch.utils.data.DataLoader(
    dataset=testset, batch_size=batch_size, shuffle=True, num_workers=1
)

x_train = trainset.data.numpy().astype(dtype="float32") / 255.0
x_test = testset.data.numpy().astype(dtype="float32") / 255.0
y_train = trainset.targets.numpy()
y_test = testset.targets.numpy()

class_per_label_size = 100
sampled_class_ids = np.concatenate(
    [
        [
            np.random.choice(
                a=np.arange(start=0, stop=y_train.shape[0], step=1)[
                    y_train == class_label
                ],
                size=class_per_label_size,
                replace=False,
            )
        ]
        for class_label in np.unique(ar=y_train)
    ],
    axis=0,
)
sampled_x_train = np.concatenate(
    [x_train[these_class_ids, :, :] for these_class_ids in sampled_class_ids], axis=0
)
sampled_y_train = np.concatenate(
    [y_train[these_class_ids] for these_class_ids in sampled_class_ids], axis=0
)<jupyter_output><empty_output><jupyter_text>MNIST: No AutoencodersWe run the MNIST dataset without using an autoencoder. The 2 dimensional tensor of pixel intensities per image for MNIST images are of dimension $\mathbb{R}^{28 \times 28}$. We reshape them as a 1 dimensional tensor of dimension $\mathbb{R}^{784}$ per image. Therefore we have 784, i.e., $28 \times 28 = 784$, features for this supervised learning task per image. ResultsThe accuracy score for the MNIST classification task without autoencoders: 87.33%.<jupyter_code>pipe_base = Pipeline(
    steps=[
        ("scaler", MinMaxScaler(feature_range=(0.0, 1.0))),
        (
            "classifier",
            linear_model.LogisticRegression(max_iter=10000, random_state=666),
        ),
    ]
)

pipe_base.fit(
    X=sampled_x_train.reshape(
        sampled_x_train.shape[0], sampled_x_train.shape[1] * sampled_x_train.shape[2]
    ),
    y=sampled_y_train,
)

acc_base = pipe_base.score(
    X=x_test.reshape(x_test.shape[0], x_test.shape[1] * x_test.shape[2]), y=y_test
)

print(
    f"The accuracy score for the MNIST classification task without autoencoders: {acc_base * 100}%"
)
tsne = TSNE(n_components=2)
test_encoded_df = pd.DataFrame(
    data=tsne.fit_transform(
        X=x_test.reshape(x_test.shape[0], x_test.shape[1] * x_test.shape[2])
    ),
    columns=["Encoding_1", "Encoding_2"],
)
test_encoded_df["Target"] = y_test

cluster_plot = (
    plotnine.ggplot(data=test_encoded_df)
    + plotnine.geom_point(
        mapping=plotnine.aes(x="Encoding_1", y="Encoding_2", fill="factor(Target)"),
        size=2,
        color="black",
    )
    + plotnine.xlab(xlab="Encoding dimension 1")
    + plotnine.ylab(ylab="Encoding dimension 2")
    + plotnine.ggtitle(title="No autoencoder with 2-dimensional encoding")
    + plotnine.theme_matplotlib()
)

cluster_plot.save(
    filename="/home/hamaad/Projects/autoencoders_pytorch/no_autoencoders_plot.png",
    height=6,
    width=8,
)<jupyter_output><empty_output><jupyter_text>![](no_autoencoders_plot.png)  MNIST: PCAWe use a PCA filter that picks the number of components that explain $99\%$ of the variation. ResultsThe accuracy score for the MNIST classification task with PCA: 78.68%.<jupyter_code>pca = PCA(n_components=0.99)
pca.fit(X=x_train.reshape(x_train.shape[0], x_train.shape[1] * x_train.shape[2]))

pipe_pca = Pipeline(
    steps=[
        ("scaler", MinMaxScaler(feature_range=(0.0, 1.0))),
        (
            "classifier",
            linear_model.LogisticRegression(max_iter=10000, random_state=666),
        ),
    ]
)

pipe_pca.fit(
    X=pca.transform(
        X=sampled_x_train.reshape(
            sampled_x_train.shape[0],
            sampled_x_train.shape[1] * sampled_x_train.shape[2],
        )
    ),
    y=sampled_y_train,
)

acc_pca = pipe_pca.score(
    X=pca.transform(
        X=x_test.reshape(x_test.shape[0], x_test.shape[1] * x_test.shape[2])
    ),
    y=y_test,
)

print(
    f"The accuracy score for the MNIST classification task with PCA: {acc_pca * 100}%"
)
tsne = TSNE(n_components=2)
test_encoded_df = pd.DataFrame(
    data=tsne.fit_transform(
        X=pca.transform(
            x_test.reshape(x_test.shape[0], x_test.shape[1] * x_test.shape[2])
        )
    ),
    columns=["Encoding_1", "Encoding_2"],
)
test_encoded_df["Target"] = y_test

cluster_plot = (
    plotnine.ggplot(data=test_encoded_df)
    + plotnine.geom_point(
        mapping=plotnine.aes(x="Encoding_1", y="Encoding_2", fill="factor(Target)"),
        size=2,
        color="black",
    )
    + plotnine.xlab(xlab="Encoding dimension 1")
    + plotnine.ylab(ylab="Encoding dimension 2")
    + plotnine.ggtitle(title="PCA with 2-dimensional encoding")
    + plotnine.theme_matplotlib()
)

cluster_plot.save(
    filename="/home/hamaad/Projects/autoencoders_pytorch/pca_plot.png",
    height=6,
    width=8,
)<jupyter_output><empty_output><jupyter_text>![](pca_plot.png)  MNIST: Vanilla AutoencodersAn autoencoder is an unsupervised learning technique where the objective is to learn a set of features that can be used to reconstruct the input data.Our input data is $X \in \mathbb{R}^{N \times 784}$. An encoder function $E$ maps this to a set of $K$ features such that $E: \mathbb{R}^{N \times 784} \rightarrow \mathbb{R}^{N \times K}$. A decoder function $D$ uses the set of $K$ features to reconstruct the input data such that $D: \mathbb{R}^{N \times K} \rightarrow \mathbb{R}^{N \times 784}$. \begin{align*}&X \in \mathbb{R}^{N \times 784} \\&E: \mathbb{R}^{N \times 784} \rightarrow \mathbb{R}^{N \times K} \\&D: \mathbb{R}^{N \times K} \rightarrow \mathbb{R}^{N \times 784}\end{align*}    Lets denote the reconstructed data as $\tilde{X} = D(E(X))$. The goal is to learn the encoding and decoding functions such that we minimize the difference between the input data and the reconstructed data. An example for an objective function for this task can be the Mean Squared Error (MSE) such that $\frac{1}{N}||\tilde{X} - X||^{2}_{2}$.     We learn the encoding and decoding functions by minimizing the MSE using the parameters that define the encoding and decoding functions: The gradient of the MSE with respect to the parameters are calculated using the chain rule, i.e., backpropagation, and used to update the parameters via an optimization algorithm such as Stochastic Gradient Descent (SGD). Lets assume we have a single layer autoencoder using the Exponential Linear Unit (ELU) activation function, batch normalization, dropout and the Adaptive Moment (Adam) optimization algorithm. $B$ is the batch size, $K$ is the number of features.* **Exponential Linear Unit:** The activation function is smooth everywhere and avoids the vanishing gradient problem as the output takes on negative values when the input is negative. $\alpha$ is taken to be $1.0$.\begin{align*}H_{\alpha}(z) &= \begin{cases}&\alpha\left(\exp(z) - 1\right) \quad \text{if} \quad z < 0 \\&z \quad \text{if} \quad z \geq 0\end{cases} \\\frac{dH_{\alpha}(z)}{dz} &= \begin{cases}&\alpha\left(\exp(z)\right) \quad \text{if} \quad z < 0 \\&1 \quad \text{if} \quad z \geq 0\end{cases} \end{align*}* **Batch Normalization:** The idea is to transform the inputs into a hidden layer's activation functions. We standardize or normalize first using the mean and variance parameters on a per feature basis and then learn a set of scaling and shifting parameters on a per feature basis that transforms the data. The following equations describe this layer succintly: The parameters we learn in this layer are $\left(\mu_{j}, \sigma_{j}^2, \beta_{j}, \gamma_{j}\right) \quad \forall j \in \{1, \dots, K\}$.\begin{align*}\mu_{j} &= \frac{1}{B} \sum_{i=1}^{B} X_{i,j} \quad &\forall j \in \{1, \dots, K\} \\\sigma_{j}^2 &= \frac{1}{B} \sum_{i=1}^{B} \left(X_{i,j} - \mu_{j}\right)^2 \quad &\forall j \in \{1, \dots, K\} \\\hat{X}_{:,j} &= \frac{X_{:,j} - \mu_{j}}{\sqrt{\sigma_{j}^2 + \epsilon}} \quad &\forall j \in \{1, \dots, K\} \\Z_{:,j} &= \gamma_{j}\hat{X}_{:,j} + \beta_{j} \quad &\forall j \in \{1, \dots, K\}\end{align*}* **Dropout:** This regularization technique simply drops the outputs from input and hidden units with a certain probability say $50\%$. * **Adam Optimization Algorithm:** This adaptive algorithm combines ideas from the Momentum and RMSProp optimization algorithms. The goal is to have some memory of past gradients which can guide future parameters updates. The following equations for the algorithm succintly describe this method assuming $\theta$ is our set of parameters to be learnt and $\eta$ is the learning rate.\begin{align*}m &\leftarrow \beta_{1}m + \left[\left(1 - \beta_{1}\right)\left(\nabla_{\theta}\text{MSE}\right)\right] \\s &\leftarrow \beta_{2}s + \left[\left(1 - \beta_{2}\right)\left(\nabla_{\theta}\text{MSE} \otimes \nabla_{\theta}\text{MSE} \right)\right] \\\theta &\leftarrow \theta - \eta m \oslash \sqrt{s + \epsilon}\end{align*} ResultsThe accuracy score for the MNIST classification task with an autoencoder: 86.67%.<jupyter_code>@typeguard.typechecked
class VanillaAutoencoder(torch.nn.Module):
    def __init__(self, num_feat: int, num_hidden: int, num_encoding: int, dropout_rate: float) -> None:
        super().__init__()

        args, _, _, values = inspect.getargvalues(frame=inspect.currentframe())
        values.pop("self")

        for arg, val in values.items():
            setattr(self, arg, val)

        self.encoder = torch.nn.Sequential(
            torch.nn.BatchNorm1d(num_features=self.num_feat, track_running_stats=True),
            torch.nn.Linear(
                in_features=self.num_feat, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_encoding, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
        )

        self.decoder = torch.nn.Sequential(
            torch.nn.BatchNorm1d(
                num_features=self.num_encoding, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_encoding, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_feat, bias=True
            ),
            torch.nn.Sigmoid(),
        )

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        x = self.encoder(input=x)
        return self.decoder(input=x)
autoencoder = VanillaAutoencoder(
    num_feat=x_train.shape[1] * x_train.shape[2],
    num_hidden=400,
    num_encoding=200,
    dropout_rate=0.5,
).to(device=device)

optimizer = torch.optim.Adam(
    lr=1e-3, amsgrad=True, params=autoencoder.parameters(), eps=1e-7
)
num_epochs = 50

for epoch in range(num_epochs):
    epoch_train_loss = []
    epoch_test_loss = []
    for X, _ in trainloader:
        X = X.view(size=(X.size(dim=0), -1)).to(device=device)

        decoding = autoencoder.forward(x=X)
        train_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)

        optimizer.zero_grad()
        train_loss.backward(retain_graph=False)
        optimizer.step()
        epoch_train_loss.append(train_loss.item())

    for X, _ in testloader:
        with torch.no_grad():
            X = X.view(size=(X.size(dim=0), -1)).to(device=device)
            decoding = autoencoder.forward(x=X)
            test_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)
            epoch_test_loss.append(test_loss.item())

    print(f"Epoch: {epoch + 1}")
    print(f"Train loss: {np.mean(a=epoch_train_loss)}")
    print(f"Test loss: {np.mean(a=epoch_test_loss)}")
pipe_autoencoder = Pipeline(
    steps=[
        ("scaler", MinMaxScaler(feature_range=(0.0, 1.0))),
        (
            "classifier",
            linear_model.LogisticRegression(max_iter=10000, random_state=666),
        ),
    ]
)

with torch.no_grad():
    pipe_autoencoder.fit(
        X=autoencoder.encoder(
            torch.FloatTensor(
                data=sampled_x_train.reshape(
                    sampled_x_train.shape[0],
                    sampled_x_train.shape[1] * sampled_x_train.shape[2],
                )
            ).to(device=device)
        )
        .cpu()
        .numpy(),
        y=sampled_y_train,
    )

    acc_autoencoder = pipe_autoencoder.score(
        X=autoencoder.encoder(
            torch.FloatTensor(data=x_test.reshape(x_test.shape[0], -1)).to(device=device)
        )
        .cpu()
        .numpy(),
        y=y_test,
    )

print(
    f"The accuracy score for the MNIST classification task with an autoencoder: {acc_autoencoder * 100}%" 
)
tsne = TSNE(n_components=2)
with torch.no_grad():
    test_encoded_df = pd.DataFrame(
        data=tsne.fit_transform(
            X=autoencoder.encoder(
                torch.FloatTensor(data=x_test.reshape(x_test.shape[0], -1)).to(device=device)
            )
            .cpu()
            .numpy()
        ),
        columns=["Encoding_1", "Encoding_2"],
    )
test_encoded_df["Target"] = y_test

cluster_plot = (
    plotnine.ggplot(data=test_encoded_df)
    + plotnine.geom_point(
        mapping=plotnine.aes(x="Encoding_1", y="Encoding_2", fill="factor(Target)"),
        size=2,
        color="black",
    )
    + plotnine.xlab(xlab="Encoding dimension 1")
    + plotnine.ylab(ylab="Encoding dimension 2")
    + plotnine.ggtitle(title="Autoencoder with 2-dimensional encoding")
    + plotnine.theme_matplotlib()
)
cluster_plot.save(
    filename="/home/hamaad/Projects/autoencoders_pytorch/autoencoders_plot.png",
    height=6,
    width=8,
)<jupyter_output><empty_output><jupyter_text>![](autoencoders_plot.png)  MNIST: Denoising AutoencodersThe idea here is to add some noise to the data and try to learn a set of robust features that can reconstruct the non-noisy data from the noisy data. The MSE objective functions is as follows, $\frac{1}{N}||D(E(X + \epsilon)) - X||^{2}_{2}$, where $\epsilon$ is some noise term.\begin{align*}&X \in \mathbb{R}^{N \times 784} \\&E: \mathbb{R}^{N \times 784} \rightarrow \mathbb{R}^{N \times K} \\&D: \mathbb{R}^{N \times K} \rightarrow \mathbb{R}^{N \times 784}\end{align*} ResultsThe accuracy score for the MNIST classification task with a denoising autoencoder: 86.62%.<jupyter_code>denoising_autoencoder = VanillaAutoencoder(
    num_feat=x_train.shape[1] * x_train.shape[2],
    num_hidden=400,
    num_encoding=200,
    dropout_rate=0.5,
).to(device=device)

optimizer = torch.optim.Adam(
    lr=1e-3, amsgrad=True, params=denoising_autoencoder.parameters(), eps=1e-7
)
num_epochs = 50

for epoch in range(num_epochs):
    epoch_train_loss = []
    epoch_test_loss = []
    for X, _ in trainloader:
        X = X.view(size=(X.size(dim=0), -1)).to(device=device)
        noise = torch.rand(size=(X.size(dim=0), X.size(dim=1))).cuda()

        decoding = denoising_autoencoder.forward(x=X + (0.1 * noise))
        train_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)

        optimizer.zero_grad()
        train_loss.backward(retain_graph=False)
        optimizer.step()
        epoch_train_loss.append(train_loss.item())

    for X, _ in testloader:
        with torch.no_grad():
            X = X.view(size=(X.size(dim=0), -1)).to(device=device)
            decoding = denoising_autoencoder.forward(x=X + (0.1 * noise))
            test_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)
            epoch_test_loss.append(test_loss.item())

    print(f"Epoch: {epoch + 1}")
    print(f"Train loss: {np.mean(a=epoch_train_loss)}")
    print(f"Test loss: {np.mean(a=epoch_test_loss)}")
pipe_denoising_autoencoder = Pipeline(
    steps=[
        ("scaler", MinMaxScaler(feature_range=(0.0, 1.0))),
        (
            "classifier",
            linear_model.LogisticRegression(max_iter=10000, random_state=666),
        ),
    ]
)

with torch.no_grad():
    pipe_denoising_autoencoder.fit(
        X=denoising_autoencoder.encoder(
            torch.FloatTensor(
                data=sampled_x_train.reshape(
                    sampled_x_train.shape[0],
                    sampled_x_train.shape[1] * sampled_x_train.shape[2],
                )
            ).to(device=device)
        )
        .cpu()
        .numpy(),
        y=sampled_y_train,
    )

    acc_denoising_autoencoder = pipe_denoising_autoencoder.score(
        X=denoising_autoencoder.encoder(
            torch.FloatTensor(data=x_test.reshape(x_test.shape[0], -1)).to(device=device)
        )
        .cpu()
        .numpy(),
        y=y_test,
    )

print(
    f"The accuracy score for the MNIST classification task with a denoising autoencoder: {acc_denoising_autoencoder * 100}"
)
tsne = TSNE(n_components=2)
with torch.no_grad():
    test_encoded_df = pd.DataFrame(
        data=tsne.fit_transform(
            X=denoising_autoencoder.encoder(
                torch.FloatTensor(data=x_test.reshape(x_test.shape[0], -1)).to(device=device)
            )
            .cpu()
            .numpy()
        ),
        columns=["Encoding_1", "Encoding_2"],
    )
test_encoded_df["Target"] = y_test

cluster_plot = (
    plotnine.ggplot(data=test_encoded_df)
    + plotnine.geom_point(
        mapping=plotnine.aes(x="Encoding_1", y="Encoding_2", fill="factor(Target)"),
        size=2,
        color="black",
    )
    + plotnine.xlab(xlab="Encoding dimension 1")
    + plotnine.ylab(ylab="Encoding dimension 2")
    + plotnine.ggtitle(title="Denoising autoencoder with 2-dimensional encoding")
    + plotnine.theme_matplotlib()
)
cluster_plot.save(
    filename="/home/hamaad/Projects/autoencoders_pytorch/denoising_autoencoders_plot.png",
    height=6,
    width=8,
)<jupyter_output><empty_output><jupyter_text>![](denoising_autoencoders_plot.png)  MNIST: 1 Dimensional Convolutional AutoencodersSo far we have used flattened or reshaped raw data. Such a 1 dimensional tensor of pixel intensities per image, $\mathbb{R}^{784}$, might not take into account useful spatial features that the 2 dimensional tensor, $\mathbb{R}^{28\times28}$, might contain. To overcome this problem, we introduce the concept of convolution filters, considering first their 1 dimensional version and then their 2 dimensional version. \begin{align*}&X \in \mathbb{R}^{N \times 28 \times 28} \\&E: \mathbb{R}^{N \times 28 \times 28} \rightarrow \mathbb{R}^{N \times K} \\&D: \mathbb{R}^{N \times K} \rightarrow \mathbb{R}^{N \times 28 \times 28}\end{align*}The ideas behind convolution filters are closely related to handcrafted feature engineering: One can view the handcrafted features as simply the result of a predefined convolution filter, i.e., a convolution filter that has not been learnt based on the raw data at hand. Suppose we have raw transactions data per some unit of analysis, i.e., mortgages, that will potentially help us in classifying a unit as either defaulted or not defaulted. We will keep this example simple by only allowing the transaction values to be either \\$100 or \\$0. The raw data per unit spans 5 time periods while the defaulted label is for the next period, i.e., period 6. Here is an example of a raw data for a particular unit:\begin{align*}x = \begin{array}{l}\text{Period 1} \\ \text{Period 2} \\ \text{Period 3} \\ \text{Period 4} \\ \text{Period 5}\end{array}    \left[    \begin{array}    {c}    \$0 \\ \$0 \\ \$100 \\ \$0 \\ \$0    \end{array}    \right]\end{align*}Suppose further that if the average transaction value is \\$20 then we will see a default in period 6 for this particular mortgage unit. Otherwise we do not see a default in period 6. The average transaction value is an example of a handcrafted feature: A predefined handcrafted feature that has not been learnt in any manner. It has been arrived at via domain knowledge of credit risk. Denote this as $\mathbf{H}(x)$.The idea of learning such a feature is an example of a 1 dimensional convolution filter. As follows:\begin{align*}\mathbf{C}(x|\alpha) = \alpha_1 x_1 + \alpha_2 x_2 + \alpha_3 x_3 + \alpha_4 x_4 + \alpha_5 x_5\end{align*}Assuming that $\mathbf{H}(x)$ is the correct representation of the raw data for this supervised learning task then the optimal set of parameters learnt via supervised learning, or perhaps unsupervised learning and then transferred to the supervised learning task, i.e., transfer learning, for $\mathbf{C}(x|\alpha)$ is as follows where $\alpha$ is $\left[0.2, 0.2, 0.2, 0.2, 0.2\right]$:\begin{align*}\mathbf{C}(x|\alpha) = 0.2 x_1 + 0.2 x_2 + 0.2 x_3 + 0.2 x_4 + 0.2 x_5\end{align*}This is a simple example however this clearly illusrates the principle behind using deep learning for automatic feature engineering or representation learning. One of the main benefits of learning such a representation in an unsupervised manner is that the same representation can then be used for multiple supervised learning tasks: Transfer learning. This is a principled manner of learning a representation from raw data.To summarize the 1 dimensional convolution filter for our simple example is defined as: \begin{align*}\mathbf{C}(x|\alpha)&= x * \alpha \\&= \sum_{t=1}^{5} x_t \alpha_t\end{align*}* $x$ is the input.* $\alpha$ is the kernel.* The output $x * \alpha$ is called a feature map and $*$ is the convolution operator or filter. This is the main difference between a vanilla neural network and a convolution neural network: We replace the matrix multiplication operator by the convolution operator.* Depending on the task at hand we can have different types of convolution filters.* Kernel size can be altered. In our example the kernel size is 5.* Stride size can be altered. In our example we had no stride size however suppose that stride size was 1 and kernel size was 2, i.e., $\alpha = \left[\alpha_1, \alpha_2\right]$, then we would apply the kernel $\alpha$ at the start of the input, i.e., $\left[x_1, x_2\right] * \left[\alpha_1, \alpha_2\right]$, and move the kernel over the next area of the input, i.e., $\left[x_2, x_3\right] * \left[\alpha_1, \alpha_2\right]$, and so on and so forth until we arrive at a feature map that consists of 4 real values. This is called a valid convolution while a padded, i.e., say padded with zero values, convolution would give us a feature map that is the same size as the input, i.e., 5 real values in our example.* We can apply an activation function to the feature maps such as ELU mentioned earlier.* Finally we can summarize the information contained in feature maps by taking a maximum or average value over a defined portion of the feature map. For instance, if after using a valid convolution we arrive at a feature map of size 4 and then apply a max pooling operation with size 4 then we will be taking the maximum value of this feature map. The result is another feature map.This automates feature engineering however introduces architecture engineering where different architectures consisting of various convolution filters, activation functions, batch normalization layers, dropout layers and pooling operators can be stacked together in a pipeline in order to learn a good representation of the raw data. One usually creates an ensemble of such architectures.The goal behind convolutional autoencoders is to use convolution filters, activation functions, batch normalization layers, dropout layers and pooling operators to create an encoder function which will learn a good representation of our raw data. The decoder will also use a similar set of layers as the encoder to reconstruct the raw data with one exception: Instead of using a pooling operator it will use an upsampling operator. The basic idea behind the upsampling operator is to repeat an element a certain number of times say size 4: One can view this as the inverse operator to the pooling operator. The pooling operator is essentially a downsampling operator and the upsampling operator is simply the inverse of that in some sense. ResultsThe accuracy score for the MNIST classification task with a 1 dimensional convolutional autoencoder: 90.08%.<jupyter_code>class ConvolutionalAutoencoder(torch.nn.Module):
    def __init__(
        self,
        filters,
        kernel_size,
        strides,
        pool_size,
        num_feat,
        dropout_rate,
    ):
        super().__init__()

        args, _, _, values = inspect.getargvalues(frame=inspect.currentframe())
        values.pop("self")

        for arg, val in values.items():
            setattr(self, arg, val)

        self.encoder = torch.nn.Sequential(
            torch.nn.BatchNorm1d(num_features=self.num_feat, track_running_stats=True),
            torch.nn.Conv1d(
                in_channels=self.num_feat,
                out_channels=self.filters,
                kernel_size=self.kernel_size,
                stride=self.strides,
                dilation=1,
                bias=True,
                padding=int((self.kernel_size - 1) / 2),
            ),
            torch.nn.ConstantPad1d(padding=(1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(num_features=self.filters, track_running_stats=True),
            torch.nn.Conv1d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=self.kernel_size,
                stride=self.strides,
                dilation=1,
                bias=True,
                padding=int((self.kernel_size - 1) / 2),
            ),
            torch.nn.ConstantPad1d(padding=(1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(num_features=self.filters, track_running_stats=True),
            torch.nn.Conv1d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=self.kernel_size,
                stride=self.strides,
                dilation=1,
                bias=True,
                padding=int((self.kernel_size - 1) / 2),
            ),
            torch.nn.ConstantPad1d(padding=(1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.MaxPool1d(
                kernel_size=self.pool_size,
                stride=self.pool_size,
                padding=0,
                dilation=1,
                return_indices=False,
                ceil_mode=False,
            ),
        )

        self.decoder = torch.nn.Sequential(
            torch.nn.BatchNorm1d(num_features=self.filters, track_running_stats=True),
            torch.nn.Conv1d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=self.kernel_size,
                stride=self.strides,
                dilation=1,
                bias=True,
                padding=int((self.kernel_size - 1) / 2),
            ),
            torch.nn.ConstantPad1d(padding=(1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Upsample(scale_factor=self.pool_size),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(num_features=self.filters, track_running_stats=True),
            torch.nn.Conv1d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=self.kernel_size,
                stride=self.strides,
                dilation=1,
                bias=True,
                padding=int((self.kernel_size - 1) / 2),
            ),
            torch.nn.ConstantPad1d(padding=(1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(num_features=self.filters, track_running_stats=True),
            torch.nn.Conv1d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=self.kernel_size,
                stride=self.strides,
                dilation=1,
                bias=True,
                padding=int((self.kernel_size - 1) / 2),
            ),
            torch.nn.ConstantPad1d(padding=(1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.BatchNorm1d(num_features=self.filters, track_running_stats=True),
            torch.nn.Conv1d(
                in_channels=self.filters,
                out_channels=self.num_feat,
                kernel_size=self.kernel_size,
                stride=self.strides,
                dilation=1,
                bias=True,
                padding=int((self.kernel_size - 1) / 2),
            ),
            torch.nn.ConstantPad1d(padding=(1, 0), value=0),
            torch.nn.Sigmoid(),
        )

    def forward(self, x):
        x = self.encoder(x)
        return self.decoder(x)
convolutional_autoencoder = ConvolutionalAutoencoder(
    num_feat=x_train.shape[1],
    filters=50,
    kernel_size=8,
    strides=1,
    pool_size=4,
    dropout_rate=0.5,
).to(device)

optimizer = torch.optim.Adam(
    lr=1e-3, amsgrad=True, params=convolutional_autoencoder.parameters(), eps=1e-7
)
num_epochs = 50

for epoch in range(num_epochs):
    epoch_train_loss = []
    epoch_test_loss = []
    for X, _ in trainloader:
        X = X.view(size=(X.size(dim=0), X.size(dim=2), X.size(dim=3))).to(device)
        decoding = convolutional_autoencoder.forward(x=X)
        train_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)

        optimizer.zero_grad()
        train_loss.backward(retain_graph=False)
        optimizer.step()
        epoch_train_loss.append(train_loss.item())

    for X, _ in testloader:
        with torch.no_grad():
            X = X.view(size=(X.size(dim=0), X.size(dim=2), X.size(dim=3))).to(device)
            decoding = convolutional_autoencoder.forward(x=X)
            test_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)
            epoch_test_loss.append(test_loss.item())

    print("Epoch:", epoch + 1)
    print("Train loss:", np.mean(a=epoch_train_loss))
    print("Test loss:", np.mean(a=epoch_test_loss))
pipe_convolutional_autoencoder = Pipeline(
    steps=[
        ("scaler", MinMaxScaler(feature_range=(0.0, 1.0))),
        (
            "classifier",
            linear_model.LogisticRegression(max_iter=10000, random_state=666),
        ),
    ]
)

with torch.no_grad():
    pipe_convolutional_autoencoder.fit(
        X=torch.nn.Flatten()(
            convolutional_autoencoder.encoder(
                torch.FloatTensor(data=sampled_x_train).to(device)
            )
        )
        .cpu()
        .numpy(),
        y=sampled_y_train,
    )

    acc_convolutional_autoencoder = pipe_convolutional_autoencoder.score(
        X=torch.nn.Flatten()(
            convolutional_autoencoder.encoder(torch.FloatTensor(data=x_test).to(device))
        )
        .cpu()
        .numpy(),
        y=y_test,
    )

print(
    "The accuracy score for the MNIST classification task with a convolutional autoencoder: %.6f%%."
    % (acc_convolutional_autoencoder * 100)
)
tsne = TSNE(n_components=2)
with torch.no_grad():
    test_encoded_df = pd.DataFrame(
        data=tsne.fit_transform(
            X=torch.nn.Flatten()(
                convolutional_autoencoder.encoder(
                    torch.FloatTensor(data=x_test).to(device)
                )
            )
            .cpu()
            .numpy()
        ),
        columns=["Encoding_1", "Encoding_2"],
    )
test_encoded_df["Target"] = y_test

cluster_plot = (
    plotnine.ggplot(data=test_encoded_df)
    + plotnine.geom_point(
        mapping=plotnine.aes(x="Encoding_1", y="Encoding_2", fill="factor(Target)"),
        size=2,
        color="black",
    )
    + plotnine.xlab(xlab="Encoding dimension 1")
    + plotnine.ylab(ylab="Encoding dimension 2")
    + plotnine.ggtitle(title="Convolutional autoencoder with 2-dimensional encoding")
    + plotnine.theme_matplotlib()
)
cluster_plot.save(
    filename="/home/hamaad/Projects/autoencoders_pytorch/convolutional_autoencoders_plot.png",
    height=6,
    width=8,
)<jupyter_output><empty_output><jupyter_text>![](convolutional_autoencoders_plot.png)  MNIST: Sequence to Sequence AutoencodersGiven our mortgage default example a potentially more useful deep learning architecture might be the Recurrent Neural Network (RNN), specifically their state of the art variant the Long Short Term Memory (LSTM) network. The goal is to explicitly take into account the sequential nature of the raw data.\begin{align*}&X \in \mathbb{R}^{N \times 28 \times 28} \\&E: \mathbb{R}^{N \times 28 \times 28} \rightarrow \mathbb{R}^{N \times K} \\&D: \mathbb{R}^{N \times K} \rightarrow \mathbb{R}^{N \times 28 \times 28}\end{align*}The gradients in a RNN depend on the parameter matrices defined for the model. Simply put these parameter matrices can end up being multiplied many times over and hence cause two major problems for learning: Exploding and vanishing gradients. If the spectral radius of the parameter matrices, i.e., the maximum absolute value of the eigenvalues of a matrix, is more than 1 then gradients can become large enough, i.e., explode in value, such that learning diverges and similarly if the spectral radius is less than 1 then gradients can become small, i.e., vanish in value, such that the next best transition for the parameters cannot be reliably calculated. Appropriate calculation of the gradient is important for estimating the optimal set of parameters that define a machine learning method and the LSTM network overcomes these problems in a vanilla RNN. We now define the LSTM network for 1 time step, i.e., 1 memory cell.We calculate the value of the input gate, the value of the memory cell state at time period $t$ where $f(x)$ is some activation function and the value of the forget gate:\begin{align*}i_{t} &= \sigma(W_{i}x_{t} + U_{i}h_{t-1} + b_{i}) \\\tilde{c_{t}} &= f(W_{c}x_{t} + U_{c}h_{t-1} + b_{c}) \\f_{t} &= \sigma(W_{f}x_{t} + U_{f}h_{t-1} + b_{f})\end{align*}The forget gate controls the amount the LSTM remembers, i.e., the value of the memory cell state at time period $t-1$ where $\otimes$ is the hadamard product:\begin{align*}c_{t} = i_{t} \otimes \tilde{c_{t}} + f_{t} \otimes c_{t-1} \end{align*}With the updated state of the memory cell we calculate the value of the outputs gate and finally the output value itself:\begin{align*}o_{t} &= \sigma(W_{o}x_{t} + U_{o}h_{t-1} + b_{o}) \\h_{t} &= o_{t} \otimes f(c_{t})\end{align*}We can have a wide variety of LSTM architectures such as the convolutional LSTM where note that we replace the matrix multiplication operators in the input gate, the initial estimate $\tilde{c_{t}}$ of the memory cell state, the forget gate and the output gate by the convolution operator $*$:\begin{align*}i_{t} &= \sigma(W_{i} * x_{t} + U_{i} * h_{t-1} + b_{i}) \\\tilde{c_{t}} &= f(W_{c} * x_{t} + U_{c} * h_{t-1} + b_{c}) \\f_{t} &= \sigma(W_{f} * x_{t} + U_{f} * h_{t-1} + b_{f}) \\c_{t} &= i_{t} \otimes \tilde{c_{t}} + f_{t} \otimes c_{t-1} \\ o_{t} &= \sigma(W_{o} * x_{t} + U_{o} * h_{t-1} + b_{o}) \\h_{t} &= o_{t} \otimes f(c_{t})\end{align*}Another popular variant is the peephole LSTM where the gates are allowed to peep at the memory cell state:\begin{align*}i_{t} &= \sigma(W_{i}x_{t} + U_{i}h_{t-1} + V_{i}c_{t-1} + b_{i}) \\\tilde{c_{t}} &= f(W_{c}x_{t} + U_{c}h_{t-1} + V_{c}c_{t-1} + b_{c}) \\f_{t} &= \sigma(W_{f}x_{t} + U_{f}h_{t-1} + V_{f}c_{t-1} + b_{f}) \\c_{t} &= i_{t} \otimes \tilde{c_{t}} + f_{t} \otimes c_{t-1} \\ o_{t} &= \sigma(W_{o}x_{t} + U_{o}h_{t-1} + V_{o}c_{t} + b_{o}) \\h_{t} &= o_{t} \otimes f(c_{t})\end{align*}The goal for the sequence to sequence autoencoder is to create a representation of the raw data using a LSTM as an encoder. This representation will be a sequence of vectors say, $h_{1}, \dots, h_{T}$, learnt from a sequence of raw data vectors say, $x_{1}, \dots, x_{T}$. The final vector of the representation, $h_{T}$, is our encoded representation, also called a context vector. This context vector is repeated as many times as the length of the sequence such that it can be used as an input to a decoder which is yet another LSTM. The decoder LSTM will use this context vector to recontruct the sequence of raw data vectors, $\tilde{x_{1}}, \dots, \tilde{x_{T}}$. If the context vector is useful in the recontruction task then it can be further used for other tasks such as predicting default risk as given in our example. ResultsThe accuracy score for the MNIST classification task with a sequence to sequence autoencoder: 96.18%.<jupyter_code>class Seq2SeqAutoencoder(torch.nn.Module):
    def __init__(
        self,
        num_feat,
        num_hidden,
        num_layers,
        dropout_rate,
    ):
        super().__init__()

        args, _, _, values = inspect.getargvalues(frame=inspect.currentframe())
        values.pop("self")

        for arg, val in values.items():
            setattr(self, arg, val)

        self.encoder = torch.nn.LSTM(
            input_size=self.num_feat,
            hidden_size=self.num_hidden,
            num_layers=self.num_layers,
            batch_first=True,
            bias=True,
            dropout=self.dropout_rate,
            bidirectional=False,
        )

        self.intermediate = torch.nn.LSTM(
            input_size=self.num_hidden,
            hidden_size=self.num_hidden,
            num_layers=self.num_layers - 1,
            batch_first=True,
            bias=True,
            dropout=self.dropout_rate,
            bidirectional=False,
        )

        self.decoder = torch.nn.LSTM(
            input_size=self.num_hidden,
            hidden_size=self.num_feat,
            num_layers=1,
            batch_first=True,
            bias=True,
            dropout=0.0,
            bidirectional=False,
        )

    def forward(self, x):
        encoding, (_, _) = self.encoder(x)
        encoding = (
            encoding[:, -1, :].unsqueeze(dim=1).repeat(1, encoding.size(dim=1), 1)
        )

        intermediate, (_, _) = self.intermediate(encoding)
        decoding, (_, _) = self.decoder(intermediate)

        return decoding
seq2seq_autoencoder = Seq2SeqAutoencoder(
    num_feat=x_train.shape[1],
    num_hidden=200,
    num_layers=3,
    dropout_rate=0.0,
).to(device)

optimizer = torch.optim.Adam(
    lr=1e-3,
    amsgrad=True,
    eps=1e-7,
    params=seq2seq_autoencoder.parameters(),
)
num_epochs = 50

for epoch in range(num_epochs):
    epoch_train_loss = []
    epoch_test_loss = []
    for X, _ in trainloader:
        X = X.view(size=(X.size(dim=0), X.size(dim=2), X.size(dim=3))).to(device)
        decoding = seq2seq_autoencoder.forward(x=X)
        train_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)

        optimizer.zero_grad()
        train_loss.backward(retain_graph=False)
        optimizer.step()
        epoch_train_loss.append(train_loss.item())

    for X, _ in testloader:
        with torch.no_grad():
            X = X.view(size=(X.size(dim=0), X.size(dim=2), X.size(dim=3))).to(device)
            decoding = seq2seq_autoencoder.forward(x=X)
            test_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)
            epoch_test_loss.append(test_loss.item())

    print("Epoch:", epoch + 1)
    print("Train loss:", np.mean(a=epoch_train_loss))
    print("Test loss:", np.mean(a=epoch_test_loss))
pipe_seq2seq_autoencoder = Pipeline(
    steps=[
        ("scaler", MinMaxScaler(feature_range=(0.0, 1.0))),
        (
            "classifier",
            linear_model.LogisticRegression(max_iter=10000, random_state=666),
        ),
    ]
)

with torch.no_grad():
    pipe_seq2seq_autoencoder.fit(
        X=seq2seq_autoencoder.encoder(
            torch.FloatTensor(data=sampled_x_train).to(device)
        )[0][:, -1, :]
        .cpu()
        .numpy(),
        y=sampled_y_train,
    )

    acc_seq2seq_autoencoder = pipe_seq2seq_autoencoder.score(
        X=seq2seq_autoencoder.encoder(torch.FloatTensor(data=x_test).to(device))[0][
            :, -1, :
        ]
        .cpu()
        .numpy(),
        y=y_test,
    )

print(
    "The accuracy score for the MNIST classification task with a sequence to sequence autoencoder: %.6f%%."
    % (acc_seq2seq_autoencoder * 100)
)
tsne = TSNE(n_components=2)
with torch.no_grad():
    test_encoded_df = pd.DataFrame(
        data=tsne.fit_transform(
            X=seq2seq_autoencoder.encoder(torch.FloatTensor(data=x_test).to(device))[0][
                :, -1, :
            ]
            .cpu()
            .numpy()
        ),
        columns=["Encoding_1", "Encoding_2"],
    )
test_encoded_df["Target"] = y_test

cluster_plot = (
    plotnine.ggplot(data=test_encoded_df)
    + plotnine.geom_point(
        mapping=plotnine.aes(x="Encoding_1", y="Encoding_2", fill="factor(Target)"),
        size=2,
        color="black",
    )
    + plotnine.xlab(xlab="Encoding dimension 1")
    + plotnine.ylab(ylab="Encoding dimension 2")
    + plotnine.ggtitle(
        title="Sequence to sequence autoencoder with 2-dimensional encoding"
    )
    + plotnine.theme_matplotlib()
)
cluster_plot.save(
    filename="/home/hamaad/Projects/autoencoders_pytorch/seq2seq_autoencoders_plot.png",
    height=6,
    width=8,
)<jupyter_output><empty_output><jupyter_text>![](seq2seq_autoencoders_plot.png)  MNIST: Variational AutoencodersWe now combine Bayesian inference with deep learning by using variational inference to train a vanilla autoencoder. This moves us towards generative modelling which can have further use cases in semi-supervised learning. The other benefit of training using Bayesian inference is that we can be more robust to higher capacity deep learners, i.e., avoid overfitting. \begin{align*}&X \in \mathbb{R}^{N \times 784} \\&E: \mathbb{R}^{N \times 784} \rightarrow \mathbb{R}^{N \times K} \\&D: \mathbb{R}^{N \times K} \rightarrow \mathbb{R}^{N \times 784}\end{align*}* Assume $X$ is our raw data while $Z$ is our learnt representation. * We have a prior belief on our learnt representation: \begin{align*}p(Z)\end{align*}* The posterior distribution for our learnt representation is: \begin{align*}p(Z|X)=\frac{p(X|Z)p(Z)}{p(X)}\end{align*}* The marginal likelihood, $p(X)$, is often intractable causing the posterior distribution, $p(Z|X)$, to be intractable:\begin{align*}p(X)=\int_{Z}p(X|Z)p(Z)dZ\end{align*}* We therefore need an approximate posterior distribution via variational inference that can deal with the intractability. This additionally also provides the benefit of dealing with large scale datasets as generally Markov Chain Monte Carlo (MCMC) methods are not well suited for large scale datasets. One might also consider Laplace approximation for the approximate posterior distribution however we will stick with variational inference as it allows a richer set of approximations compared to Laplace approximation. Laplace approximation simply amounts to finding the Maximum A Posteriori (MAP) estimate to an augmented likelihood optimization, taking the negative of the inverse of the Hessian at the MAP estimate to estimate the variance-covariance matrix and finally use the variance-covariance matrix with a multivariate Gaussian distribution or some other appropriate multivariate distribution.* Assume that our approximate posterior distribution, which is also our probabilistic encoder, is given as:\begin{align*}q(Z|X)\end{align*}* Our probabilistic decoder is given by:\begin{align*}p(X|Z)\end{align*}* Given our setup above with regards to an encoder and a decoder let us now write down the optimization problem where $\theta$ are the generative model parameters while $\phi$ are the variational parameters:\begin{align*}\log{p(X)}= \underbrace{D_{KL}(q(Z|X)||p(Z|X))}_\text{Intractable as p(Z|X) is intractable} + \underbrace{\mathcal{L}(\theta, \phi|X)}_\text{Evidence Lower Bound or ELBO}\end{align*}* Note that $D_{KL}(q(Z|X)||p(Z|X))$ is non-negative therefore that makes the ELBO a lower bound on $\log{p(X)}$:\begin{align*}\log{p(X)}\geq \mathcal{L}(\theta, \phi|X) \quad \text{as} \quad D_{KL}(q(Z|X)||p(Z|X)) \geq 0\end{align*}* Therefore we can alter our optimization problem to look only at the ELBO:\begin{align*}\mathcal{L}(\theta, \phi|X) &= \mathbb{E}_{q(Z|X)}\left[\log{p(X,Z)} - \log{q(Z|X)}\right] \\&= \mathbb{E}_{q(Z|X)}\left[\underbrace{\log{p(X|Z)}}_\text{Reconstruction error} + \log{p(Z)} - \log{q(Z|X)}\right] \\&= \mathbb{E}_{q(Z|X)}\left[\underbrace{\log{p(X|Z)}}_\text{Reconstruction error} - \underbrace{D_{KL}(q(Z|X)||p(Z))}_\text{Regularization}\right] \\&= \int_{Z} \left[\log{p(X|Z)} - D_{KL}(q(Z|X)||p(Z))\right] q(Z|X) dZ\end{align*}* The above integration problem can be solved via Monte Carlo integration as $D_{KL}(q(Z|X)||p(Z))$ is not intractable. Assuming that the probabilistic encoder $q(Z|X)$ is a multivariate Gaussian with a diagonal variance-covariance matrix we use the reparameterization trick to sample from this distribution say $M$ times in order to calculate the expectation term in the ELBO optimization problem. The reparameterization trick in this particular case amounts to sampling $M$ times from the standard Gaussian distribution, multiplying the samples by $\sigma$ and adding $\mu$ to the samples.  * $\mu$ is our learnt representation used for the reconstruction of the raw data. If the learnt representation is useful it can then be used for other tasks as well.* This is a powerful manner of combining Bayesian inference with deep learning. Variational inference used in this manner can be applied to various deep learning architectures and has further links with the Generative Adversarial Network (GAN). We explore the use of adversarial learning in representation learning in another repo/paper. ResultsThe accuracy score for the MNIST classification task with a variational autoencoder: 91.80%.<jupyter_code>class VariationalAutoencoder(torch.nn.Module):
    def __init__(self, num_feat, num_hidden, num_encoding, dropout_rate):
        super().__init__()

        args, _, _, values = inspect.getargvalues(frame=inspect.currentframe())
        values.pop("self")

        for arg, val in values.items():
            setattr(self, arg, val)

        self.encoder = torch.nn.Sequential(
            torch.nn.BatchNorm1d(num_features=self.num_feat, track_running_stats=True),
            torch.nn.Linear(
                in_features=self.num_feat, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_encoding, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
        )

        self.mu = torch.nn.Linear(
            in_features=self.num_encoding, out_features=self.num_encoding, bias=True
        )
        self.log_sigma = torch.nn.Linear(
            in_features=self.num_encoding, out_features=self.num_encoding, bias=True
        )

        self.decoder = torch.nn.Sequential(
            torch.nn.BatchNorm1d(
                num_features=self.num_encoding, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_encoding, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_hidden, bias=True
            ),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Dropout(p=self.dropout_rate, inplace=False),
            torch.nn.BatchNorm1d(
                num_features=self.num_hidden, track_running_stats=True
            ),
            torch.nn.Linear(
                in_features=self.num_hidden, out_features=self.num_feat, bias=True
            ),
            torch.nn.Sigmoid(),
        )

    def forward(self, x, prior):
        encoding = self.encoder(x)
        mu = self.mu(encoding)
        log_sigma = self.log_sigma(encoding)

        posterior = mu + (torch.exp(input=log_sigma / 2) * prior)

        decoding = self.decoder(posterior)

        return decoding, mu, log_sigma

    def vae_loss(self, y_true, y_pred, mu, log_sigma, beta):
        recon = torch.nn.MSELoss(reduction="mean")(input=y_pred, target=y_true)
        kl = -0.5 * torch.sum(
            input=1 + log_sigma - torch.exp(input=log_sigma) - mu.pow(exponent=2),
            axis=-1,
        )
        return recon + (beta * torch.mean(input=kl))

    def feature_extractor(self, x):
        return self.mu(self.encoder(x))
beta = batch_size / x_train.shape[0]
variational_autoencoder = VariationalAutoencoder(
    num_feat=x_train.shape[1] * x_train.shape[2],
    num_hidden=400,
    num_encoding=200,
    dropout_rate=0.0,
).to(device)

optimizer = torch.optim.Adam(
    lr=1e-3, amsgrad=True, params=variational_autoencoder.parameters(), eps=1e-7
)
num_epochs = 50

for epoch in range(num_epochs):
    epoch_train_loss = []
    epoch_test_loss = []
    for X, _ in trainloader:
        X = X.view(size=(X.size(dim=0), -1)).to(device)
        prior = torch.randn(
            size=(X.size(dim=0), variational_autoencoder.num_encoding)
        ).to(device)

        decoding, mu, log_sigma = variational_autoencoder.forward(x=X, prior=prior)
        train_loss = variational_autoencoder.vae_loss(
            y_true=X, y_pred=decoding, mu=mu, log_sigma=log_sigma, beta=beta
        )

        optimizer.zero_grad()
        train_loss.backward(retain_graph=False)
        optimizer.step()
        epoch_train_loss.append(train_loss.item())

    for X, _ in testloader:
        with torch.no_grad():
            X = X.view(size=(X.size(dim=0), -1)).to(device)
            prior = torch.randn(
                size=(X.size(dim=0), variational_autoencoder.num_encoding)
            ).to(device)
            decoding, mu, log_sigma = variational_autoencoder.forward(x=X, prior=prior)
            test_loss = variational_autoencoder.vae_loss(
                y_true=X, y_pred=decoding, mu=mu, log_sigma=log_sigma, beta=beta
            )
            epoch_test_loss.append(test_loss.item())

    print("Epoch:", epoch + 1)
    print("Train loss:", np.mean(a=epoch_train_loss))
    print("Test loss:", np.mean(a=epoch_test_loss))
pipe_variational_autoencoder = Pipeline(
    steps=[
        ("scaler", MinMaxScaler(feature_range=(0.0, 1.0))),
        (
            "classifier",
            linear_model.LogisticRegression(max_iter=10000, random_state=666),
        ),
    ]
)

with torch.no_grad():
    pipe_variational_autoencoder.fit(
        X=variational_autoencoder.feature_extractor(
            torch.FloatTensor(
                data=sampled_x_train.reshape(
                    sampled_x_train.shape[0],
                    sampled_x_train.shape[1] * sampled_x_train.shape[2],
                )
            ).to(device)
        )
        .cpu()
        .numpy(),
        y=sampled_y_train,
    )

    acc_variational_autoencoder = pipe_variational_autoencoder.score(
        X=variational_autoencoder.feature_extractor(
            torch.FloatTensor(data=x_test.reshape(x_test.shape[0], -1)).to(device)
        )
        .cpu()
        .numpy(),
        y=y_test,
    )

print(
    "The accuracy score for the MNIST classification task with a variational autoencoder: %.6f%%."
    % (acc_variational_autoencoder * 100)
)
tsne = TSNE(n_components=2)
with torch.no_grad():
    test_encoded_df = pd.DataFrame(
        data=tsne.fit_transform(
            X=variational_autoencoder.feature_extractor(
                torch.FloatTensor(data=x_test.reshape(x_test.shape[0], -1)).to(device)
            )
            .cpu()
            .numpy()
        ),
        columns=["Encoding_1", "Encoding_2"],
    )
test_encoded_df["Target"] = y_test

cluster_plot = (
    plotnine.ggplot(data=test_encoded_df)
    + plotnine.geom_point(
        mapping=plotnine.aes(x="Encoding_1", y="Encoding_2", fill="factor(Target)"),
        size=2,
        color="black",
    )
    + plotnine.xlab(xlab="Encoding dimension 1")
    + plotnine.ylab(ylab="Encoding dimension 2")
    + plotnine.ggtitle(title="Variational autoencoder with 2-dimensional encoding")
    + plotnine.theme_matplotlib()
)
cluster_plot.save(
    filename="/home/hamaad/Projects/autoencoders_pytorch/variational_autoencoders_plot.png",
    height=6,
    width=8,
)

n = 30
data_size = 28
figure = np.zeros(shape=(data_size * n, data_size * n))

variational_autoencoder.eval()
with torch.no_grad():
    for i in range(n):
        for j in range(n):
            decoded = (
                variational_autoencoder.decoder(
                    torch.FloatTensor(
                        data=np.array(
                            object=np.random.normal(
                                loc=0.0,
                                scale=1.0,
                                size=(1, variational_autoencoder.num_encoding),
                            )
                        )
                    ).to(device)
                )
                .cpu()
                .numpy()
            )
            figure[
                i * data_size : (i + 1) * data_size,
                j * data_size : (j + 1) * data_size,
            ] = decoded.reshape(data_size, data_size)

plt.figure(figsize=(20, 20))
plt.imshow(X=figure, cmap="Greys_r")
plt.title(
    label="Variational Autoencoder (VAE) with 2-dimensional encoding\nGenerating new images"
)
plt.xlabel(xlabel="Encoding dimension 1")
plt.ylabel(ylabel="Encoding dimension 2")
plt.savefig(fname="/home/hamaad/Projects/autoencoders_pytorch/VAE_Generated_Images.png")
plt.show()<jupyter_output><empty_output><jupyter_text>![](VAE_Generated_Images.png) ![](variational_autoencoders_plot.png)  MNIST: 2 Dimensional Convolutional AutoencodersFor 2 dimensional convolution filters the idea is similar as for the 1 dimensional convolution filters. We will stick to our previously mentioned banking example to illustrate this point.\begin{align*}x = \begin{array}{l}\text{Period 1} \\ \text{Period 2} \\ \text{Period 3} \\ \text{Period 4} \\ \text{Period 5}\end{array}    \left[    \begin{array}    {ccc}    \$0 & \$0 & \$0 \\    \$0 & \$200 & \$0 \\    \$100 & \$0 & \$0 \\    \$0 & \$0 & \$300 \\    \$0 & \$0 & \$0    \end{array}    \right]\end{align*}In the 2 dimensional tensor of raw transactions data now we have 5 historical time periods, i.e., the rows, and 3 different transaction types, i.e., the columns. We will use a kernel, $\alpha \in \mathbb{R}^{2\times3}$, to extract useful features from the raw data. The choice of such a kernel means that we are interested in finding a feature map across all 3 transaction types and 2 historical time periods. We will use a stride length of 1 and a valid convolution to extract features over different patches of the raw data. The following will illustrate this point where $x_{\text{patch}} \subset x$:\begin{align*}\alpha &=    \left[    \begin{array}    {ccc}    \alpha_{1,1} & \alpha_{1,2} & \alpha_{1,3} \\    \alpha_{2,1} & \alpha_{2,2} & \alpha_{2,3}    \end{array}    \right] \\x_{\text{patch}} &=     \left[    \begin{array}    {ccc}    \$0 & \$0 & \$0 \\    \$0 & \$200 & \$0    \end{array}    \right] \\\mathbf{C}(x=x_{\text{patch}}|\alpha) &= x * \alpha \\&= \sum_{t=1}^{2} \sum_{k=1}^{3} x_{t,k} \alpha_{t,k}\end{align*}The principles and ideas apply to 2 dimensional convolution filters as they do for their 1 dimensional counterparts there we will not repeat them here.\begin{align*}&X \in \mathbb{R}^{N \times 28 \times 28} \\&E: \mathbb{R}^{N \times 28 \times 28} \rightarrow \mathbb{R}^{N \times K} \\&D: \mathbb{R}^{N \times K} \rightarrow \mathbb{R}^{N \times 28 \times 28}\end{align*} ResultsThe accuracy score for the MNIST classification task with a 2 dimensional convolutional autoencoder: 93.49%.<jupyter_code>class Convolutional2DAutoencoder(torch.nn.Module):
    def __init__(
        self,
        filters,
        kernel_size,
        strides,
        pool_size,
        num_channels,
        dropout_rate,
    ):
        super().__init__()

        args, _, _, values = inspect.getargvalues(frame=inspect.currentframe())
        values.pop("self")

        for arg, val in values.items():
            setattr(self, arg, val)

        self.encoder = torch.nn.Sequential(
            torch.nn.BatchNorm2d(
                num_features=self.num_channels, track_running_stats=True
            ),
            torch.nn.Conv2d(
                in_channels=self.num_channels,
                out_channels=self.filters,
                kernel_size=(self.kernel_size, self.kernel_size),
                stride=(self.strides, self.strides),
                dilation=(1, 1),
                bias=True,
                padding=((self.kernel_size - 1) // 2, (self.kernel_size - 1) // 2),
            ),
            torch.nn.ConstantPad2d(padding=(1, 0, 1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.BatchNorm2d(num_features=self.filters, track_running_stats=True),
            torch.nn.Dropout2d(p=self.dropout_rate, inplace=False),
            torch.nn.Conv2d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=(self.kernel_size, self.kernel_size),
                stride=(self.strides, self.strides),
                dilation=(1, 1),
                bias=True,
                padding=((self.kernel_size - 1) // 2, (self.kernel_size - 1) // 2),
            ),
            torch.nn.ConstantPad2d(padding=(1, 0, 1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.BatchNorm2d(num_features=self.filters, track_running_stats=True),
            torch.nn.Dropout2d(p=self.dropout_rate, inplace=False),
            torch.nn.Conv2d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=(self.kernel_size, self.kernel_size),
                stride=(self.strides, self.strides),
                dilation=(1, 1),
                bias=True,
                padding=((self.kernel_size - 1) // 2, (self.kernel_size - 1) // 2),
            ),
            torch.nn.ConstantPad2d(padding=(1, 0, 1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.MaxPool2d(
                kernel_size=(self.pool_size, self.pool_size),
                stride=(self.pool_size, self.pool_size),
                padding=(0, 0),
                dilation=(1, 1),
                return_indices=False,
                ceil_mode=False,
            ),
        )

        self.decoder = torch.nn.Sequential(
            torch.nn.BatchNorm2d(num_features=self.filters, track_running_stats=True),
            torch.nn.Conv2d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=(self.kernel_size, self.kernel_size),
                stride=(self.strides, self.strides),
                dilation=(1, 1),
                bias=True,
                padding=((self.kernel_size - 1) // 2, (self.kernel_size - 1) // 2),
            ),
            torch.nn.ConstantPad2d(padding=(1, 0, 1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.Upsample(scale_factor=self.pool_size),
            torch.nn.BatchNorm2d(num_features=self.filters, track_running_stats=True),
            torch.nn.Conv2d(
                in_channels=self.filters,
                out_channels=self.filters,
                kernel_size=(self.kernel_size, self.kernel_size),
                stride=(self.strides, self.strides),
                dilation=(1, 1),
                bias=True,
                padding=((self.kernel_size - 1) // 2, (self.kernel_size - 1) // 2),
            ),
            torch.nn.ConstantPad2d(padding=(1, 0, 1, 0), value=0),
            torch.nn.ELU(alpha=1.0, inplace=False),
            torch.nn.BatchNorm2d(num_features=self.filters, track_running_stats=True),
            torch.nn.Dropout2d(p=self.dropout_rate, inplace=False),
            torch.nn.Conv2d(
                in_channels=self.filters,
                out_channels=self.num_channels,
                kernel_size=(self.kernel_size, self.kernel_size),
                stride=(self.strides, self.strides),
                dilation=(1, 1),
                bias=True,
                padding=((self.kernel_size - 1) // 2, (self.kernel_size - 1) // 2),
            ),
            torch.nn.ConstantPad2d(padding=(1, 0, 1, 0), value=0),
            torch.nn.Sigmoid(),
        )

    def forward(self, x):
        encoding = self.encoder(x)

        return self.decoder(encoding)
convolutional2D_autoencoder = Convolutional2DAutoencoder(
    filters=50,
    kernel_size=8,
    strides=1,
    pool_size=4,
    num_channels=1,
    dropout_rate=0.5,
).to(device)

optimizer = torch.optim.Adam(
    lr=1e-3, amsgrad=True, params=convolutional2D_autoencoder.parameters(), eps=1e-7
)
num_epochs = 50

for epoch in range(num_epochs):
    epoch_train_loss = []
    epoch_test_loss = []
    for X, _ in trainloader:
        X = X.to(device)
        decoding = convolutional2D_autoencoder.forward(x=X)
        train_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)

        optimizer.zero_grad()
        train_loss.backward(retain_graph=False)
        optimizer.step()
        epoch_train_loss.append(train_loss.item())

    for X, _ in testloader:
        with torch.no_grad():
            X = X.to(device)
            decoding = convolutional2D_autoencoder.forward(x=X)
            test_loss = torch.nn.MSELoss(reduction="mean")(input=decoding, target=X)
            epoch_test_loss.append(test_loss.item())

    print("Epoch:", epoch + 1)
    print("Train loss:", np.mean(a=epoch_train_loss))
    print("Test loss:", np.mean(a=epoch_test_loss))
pipe_convolutional2D_autoencoder = Pipeline(
    steps=[
        ("scaler", MinMaxScaler(feature_range=(0.0, 1.0))),
        (
            "classifier",
            linear_model.LogisticRegression(max_iter=10000, random_state=666),
        ),
    ]
)

with torch.no_grad():
    pipe_convolutional2D_autoencoder.fit(
        X=torch.nn.Flatten()(
            convolutional2D_autoencoder.encoder(
                torch.FloatTensor(data=sampled_x_train).unsqueeze(dim=1).to(device)
            )
        )
        .cpu()
        .numpy(),
        y=sampled_y_train,
    )

    acc_convolutional2D_autoencoder = pipe_convolutional2D_autoencoder.score(
        X=torch.nn.Flatten()(
            convolutional2D_autoencoder.encoder(
                torch.FloatTensor(data=x_test).unsqueeze(dim=1).to(device)
            )
        )
        .cpu()
        .numpy(),
        y=y_test,
    )

print(
    "The accuracy score for the MNIST classification task with a 2 dimensional convolutional autoencoder: %.6f%%."
    % (acc_convolutional2D_autoencoder * 100)
)
tsne = TSNE(n_components=2)
with torch.no_grad():
    test_encoded_df = pd.DataFrame(
        data=tsne.fit_transform(
            X=torch.nn.Flatten()(
                convolutional2D_autoencoder.encoder(
                    torch.FloatTensor(data=x_test).unsqueeze(dim=1).to(device)
                )
            )
            .cpu()
            .numpy()
        ),
        columns=["Encoding_1", "Encoding_2"],
    )
test_encoded_df["Target"] = y_test

cluster_plot = (
    plotnine.ggplot(data=test_encoded_df)
    + plotnine.geom_point(
        mapping=plotnine.aes(x="Encoding_1", y="Encoding_2", fill="factor(Target)"),
        size=2,
        color="black",
    )
    + plotnine.xlab(xlab="Encoding dimension 1")
    + plotnine.ylab(ylab="Encoding dimension 2")
    + plotnine.ggtitle(title="2D convolutional autoencoder with 2-dimensional encoding")
    + plotnine.theme_matplotlib()
)
cluster_plot.save(
    filename="/home/hamaad/Projects/autoencoders_pytorch/convolutional2D_autoencoders_plot.png",
    height=6,
    width=8,
)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12014282410084422, 'fraction_numerical': 0.015279865849427548, 'mean_word_length': 2.8438073012168696, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 41, 'https://': 0, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '373004', 'license': 'Apache-2.0', 'n_tokens_mistral': 22458, 'n_tokens_neox': 21501, 'n_words': 5264, 'path': 'autoencoders.ipynb', 'repo_name': 'hamaadshah/autoencoders_pytorch', 'size': '69195'}