Autoencoder for dimensionality reduction of binary dataset for clustering

Given a binary dataset (derived from yes/no questionnaire responses) aimed to use for subsequent unsupervised cluster analysis, with significant multicolinearity and a total of 31 features by ~50,000 observations (subjects), it appeared sensible to reduce the dimensionality of the input data before performing cluster analysis. I attempted using an autoencoder for this, but surprisingly, the clusters (derived through k-medoids, chosen due to the nominal nature of the underlying data and a greater stability in relation to outliers/noise compared to e.g., k-means) were actually more distinct and clearly distinguished when using MCA, with a clear maximum Silhouette coefficient at k = 5.

Given that MCA with the first 5 PCs (explaining just ~75% of the variance, chosen through a scree plot) was used before I attempted the autoencoder way, it surprises me that an autoencoder did a worse job at extracting meaningful features at the same bottleneck dimension. The problem with the current autoencoder, appears to be that the data in the bottleneck layer, which is used in the clustering, is distorted...

Below is the code I used to construct the autoencoder. Can it be so that the hyper-parameters are off, or some details of the overall architecture? Random search of specific numbers of the number of layers, learning rate, batch size, dimensions in the layers etc. have not yielded anything substantial. Loss is similar between train and validation dataset, and levels out at around 0.15 after ~40 epochs.

I've also tried to identify studies where such data has been used, but not found anything useful.

import numpy as np
import tensorflow as tf
from tensorflow.keras.layers import Input, Dense
from tensorflow.keras.models import Model
from tensorflow.keras.optimizers import Adam
import matplotlib.pyplot as plt
input_dim = 31
layer_1_dim = 18
layer_2_dim = 10
bottleneck_dim = 5
learning_rate = 0.001
epochs = 100
batch_size = 300
# split data into training and validation
training_n = int(data.shape[0] * 0.8)
train_data = data[:training_n, :]
val_data = data[training_n:, :]
# define autoencoder initializer
initializer = tf.keras.initializers.GlorotUniform()
# autoencoder layers
input_layer = Input(shape=(input_dim,))
layer = Dense(layer_1_dim, activation='relu')(input_layer)
layer = Dense(layer_2_dim, activation='relu', kernel_initializer=initializer)(layer)
layer = Dense(bottleneck_dim, activation='relu', kernel_initializer=initializer, name="bottleneck-output")(layer)
layer = Dense(layer_2_dim, activation='relu', kernel_initializer=initializer)(layer) 
layer = Dense(layer_1_dim, activation='relu', kernel_initializer=initializer)(layer)
output_layer = Dense(input_dim, activation='sigmoid', kernel_initializer=initializer)(layer)
# define and compile autoencoder model
autoencoder = Model(inputs=input_layer, outputs=output_layer)
optimizer = Adam(learning_rate=learning_rate)
autoencoder.compile(optimizer=optimizer, loss='binary_crossentropy')
# train the autoencoder model
history = autoencoder.fit(train_data, train_data, epochs=epochs, batch_size=batch_size, validation_data=(val_data, val_data))
# get bottleneck output
bottleneck_autoencoder = Model(inputs=autoencoder.input, outputs=autoencoder.get_layer('bottleneck-output').output)
bottleneck_output = bottleneck_autoencoder.predict(data)
# plot loss in traning and validation set
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.title('Autoencoder loss (binary cross-entropy)')
plt.ylabel('Loss')
plt.xlabel('Epoch')
plt.legend(['Train', 'Validation'], loc='upper right')
plt.savefig('output/embedding.png')

So the problem with multicollinearity is that, on top of not allowing for meaningful inference, it hamper your convergence rate.
This implies that, everything equal, if your data is highly multicollinear you need significantly more data to achieve convergence.

According to what your objective is (pure prediction vs. inference) you might want to either orthogonalize (see here and here) your input and then feed the residuals to your autoencoder or you might want to choose a more sensible activation function or go with a sparse autoencoder.

That said, I would suggest to try the last two steps first if you are aiming to have good predictions and do not necessarily care about inference:

1. change your activation function in the bottleneck (use a Sigmoid?).
2. use a sparse autoencoder see here, here and here

Since the main problem here is multicollinearity, you can try using a dimensionality reduction method that factors in non-linear relationships within the data. One common example is kernelPCA
Then you can explore different kernel functions, I've usually had interesting results with just linear kernel.