Restricted Boltzmann Machine: W matrix visualization results after training MNIST images and Pseudo-log-likelihood

Question

I am implementing RBM from scratch using Tensorflow and after training my RBM on the MNIST dataset for 200 epochs using Persistent CD with two steps of contrastive divergence, I learn the weights W and hidden_bias and visible_bias for 500 hidden units. I also use a learning rate of 0.1 at the beginning, but after 100 epochs I decrease the learning rate to 0.01. I tried different batch sizes between 20 to 100, and visualized weights W matrix, as well as keeping track of Mean Square Error (MSE) and Pseudo-log-likelihood values defined according to Sklearn and LISA lab, both are decreasing overall with some fluctuations as anticipated.

I have two questions:

1. I am concerned that the visualization of my W matrix after 200 epochs does not look like pen strokes, but for the Lisa lab, after 200 epochs with a persistent CD with k=15 I get W visualizations that look like pen strokes (first figure here). Also tried Pydeep and their final W matrix using 200 epochs and a one-step CD looks like pen stoke as well. I have attached my W visualization, some of them look like digits, I am not sure if I am doing anything wrong.

2. My Pseudo-log-likelihood values start with ~-400 similar to the start value for Lisa's lab and Sklearn using their own defined functions to calculate it. Eventually, my Pseudo-log-likelihood decreases to values close to -1, after 200 epochs. But I see that for Lisa'lab it ends up close to -60, using their defined function. For Sklearn, it decreases to -6 so that value is closer to what I have. I am concerned if I am doing it incorrectly because there is a huge gap between my Pseudo-log-likelihood and Lisa's lab after 200 epochs. I already compared my functions with theirs, but I cannot catch any bugs in my code.

My score_sample(v) function for calculating Pseudo-log-likelihood (with reference):

def free_energy(v):
    with tf.name_scope('free_energy'):
        T1 = -tf.einsum('ij,j->i', v, b_v)
        T2 = -tf.reduce_sum(tf.nn.softplus(_propup(v) + b_h), axis=1)
        fe = tf.reduce_mean(T1 + T2, axis=0)
    return fe

def _propup(v):
    with tf.name_scope('prop_up'):
        t = tf.matmul(v, W)
    return t

def _propdown(h):
    with tf.name_scope('prop_down'):
        t = tf.matmul(a=h, b=W, transpose_b=True)
    return t


def log_logistic(X, out=None):
    return tf.math.log(tf.math.sigmoid(X))

def score_samples(v):
    # Randomly corrupt one feature in each sample in v.
    ind = (np.arange(v.shape[0]), np.random.randint(0, v.shape[1], v.shape[0]))
    v_ = v.numpy()
    v_[ind] = 1 - v_[ind]
    fe = free_energy(v)
    fe_ = free_energy(v_)
    return v.shape[1] * log_logistic(fe_ - fe)

I appreciate any insights on if what I have done is correct or how I can improve it?

I have already read all the posts on StackExchange about this issue and have had no luck in finding any answers to this specific question. Thank you in advance!