Comparizon of NumPy, Cudamat and CudaLearn proposed RBM code

cudalearn.cs

    public static class GpuRbmExample
    {
        public static void Main(string[] args)
        {
            CudaLearnModule.Initialize();

            var database = new MnistDatabase();

            // Training parameters

            float epsilon = 0.01f;
            float momentum = 0.9f;

            int num_epochs = 10;
            int batch_size = 64;

            int num_batches = database.Samples / batch_size;

            // Training parameters

            int num_vis = database.SampleSize;
            int num_hid = 1024;

            // Initialize Weights
            var w_vh = 0.1f * GpuMatrix<float>.Normal(num_vis, num_hid);
            var w_v = GpuMatrix<float>.Zeroes(num_vis, 1);
            var w_h = GpuMatrix<float>.Zeroes(num_hid, 1);

            // Initialize Weights updates
            var wu_vh = GpuMatrix<float>.Zeroes(num_vis, num_hid);
            var wu_v = GpuMatrix<float>.Zeroes(num_vis, 1);
            var wu_h = GpuMatrix<float>.Zeroes(num_hid, 1);

            var watch = new Stopwatch();

            var error = new float[num_epochs];
            for (int epoch = 0; epoch < num_epochs; epoch++)
            {
                Console.WriteLine(string.Format("Epoch {0}", epoch + 1));

                foreach (var v_true in database.Batches(batch_size))
                {
                    GpuMatrix<float> v = v_true;

                    // Apply momentum
                    wu_vh *= momentum;
                    wu_v *= momentum;
                    wu_h *= momentum;

                    // Positive phase
                    var h = 1.0f / (1 + Functions.Exp(-(Functions.Dot(w_vh.Transpose(), v) + w_h)));

                    wu_vh += Functions.Dot(v, h.Transpose());
                    wu_v += Functions.Sum(v, Axis.Columns);
                    wu_h += Functions.Sum(h, Axis.Columns);

                    // Sample hiddens
                    h = 1.0f * (h > GpuMatrix<float>.Uniform(num_hid, batch_size));

                    // Negative phase
                    v = 1.0f / (1 + Functions.Exp(-(Functions.Dot(w_vh, h) + w_v)));
                    h = 1.0f / (1 + Functions.Exp(-(Functions.Dot(w_vh.Transpose(), v) + w_h)));

                    wu_vh -= Functions.Dot(v, h.Transpose());
                    wu_v -= Functions.Sum(v, Axis.Columns);
                    wu_h -= Functions.Sum(h, Axis.Columns);

                    // Update weights
                    w_vh += epsilon / batch_size * wu_vh;
                    w_v += epsilon / batch_size * wu_v;
                    w_h += epsilon / batch_size * wu_h;

                    error[epoch] = Functions.Mean((v - v_true)^2, Axis.None);
                }
            }

            Console.WriteLine(string.Format("Mean squared error: {0}", Functions.Mean(error.AsMatrix(), Axis.None)));
            Console.WriteLine(string.Format("Time: {0}", watch.Elapsed));

            CudaLearnModule.Release();
        }
    }

cudamat.py

import time
import numpy as np
import cudamat as cm
import util

# initialize CUDA
cm.cublas_init()
cm.CUDAMatrix.init_random(1)

# load data
util.load('mnist.dat', globals())
dev_dat = cm.CUDAMatrix(cm.reformat(dat/255.))

# training parameters
epsilon = 0.1
momentum = 0.9

num_epochs = 30
batch_size = 128
num_batches = dat.shape[1]/batch_size

# model parameters
num_vis = dat.shape[0]
num_hid = 4096

# initialize weights
w_vh = cm.CUDAMatrix(0.1 * np.random.randn(num_vis, num_hid))
w_v = cm.CUDAMatrix(np.zeros((num_vis, 1)))
w_h = cm.CUDAMatrix(-4.*np.ones((num_hid, 1)))

# initialize weight updates
wu_vh = cm.CUDAMatrix(np.zeros((num_vis, num_hid)))
wu_v = cm.CUDAMatrix(np.zeros((num_vis, 1)))
wu_h = cm.CUDAMatrix(np.zeros((num_hid, 1)))

# initialize temporary storage
v = cm.empty((num_vis, batch_size))
h = cm.empty((num_hid, batch_size))
r = cm.empty((num_hid, batch_size))

start_time = time.time()
for epoch in range(num_epochs):
    print "Epoch " + str(epoch + 1)
    err = []

    for batch in range(num_batches):
        # get current minibatch
        v_true = dev_dat.slice(batch*batch_size,(batch + 1)*batch_size)
        v.assign(v_true)

        # apply momentum
        wu_vh.mult(momentum)
        wu_v.mult(momentum)
        wu_h.mult(momentum)

        # positive phase
        cm.dot(w_vh.T, v, target = h)
        h.add_col_vec(w_h)
        h.apply_sigmoid()

        wu_vh.add_dot(v, h.T)
        wu_v.add_sums(v, axis = 1)
        wu_h.add_sums(h, axis = 1)

        # sample hiddens
        r.fill_with_rand()
        r.less_than(h, target = h)

        # negative phase
        cm.dot(w_vh, h, target = v)
        v.add_col_vec(w_v)
        v.apply_sigmoid()

        cm.dot(w_vh.T, v, target = h)
        h.add_col_vec(w_h)
        h.apply_sigmoid()

        wu_vh.subtract_dot(v, h.T)
        wu_v.add_sums(v, axis = 1, mult = -1.)
        wu_h.add_sums(h, axis = 1, mult = -1.)

        # update weights
        w_vh.add_mult(wu_vh, epsilon/batch_size)
        w_v.add_mult(wu_v, epsilon/batch_size)
        w_h.add_mult(wu_h, epsilon/batch_size)

        # calculate reconstruction error
        v.subtract(v_true)
        err.append(v.euclid_norm()**2/(num_vis*batch_size))

    print "Mean squared error: " + str(np.mean(err))
    print "Time: " + str(time.time() - start_time)

w_vh.copy_to_host()
util.save('weights.dat', 'w_vh', {'w_vh': w_vh.numpy_array})

cm.cublas_shutdown()

numpy.py

import time
import numpy as np
import util

# load data
util.load('mnist.dat', globals())
dat = dat/255.

# training parameters
epsilon = 0.01
momentum = 0.9

num_epochs = 10
batch_size = 64
num_batches = dat.shape[1]/batch_size

# model parameters
num_vis = dat.shape[0]
num_hid = 1024

# initialize weights
w_vh = 0.1 * np.random.randn(num_vis, num_hid)
w_v = np.zeros((num_vis, 1))
w_h = np.zeros((num_hid, 1))

# initialize weight updates
wu_vh = np.zeros((num_vis, num_hid))
wu_v = np.zeros((num_vis, 1))
wu_h = np.zeros((num_hid, 1))

start_time = time.time()
for epoch in range(num_epochs):
    print "Epoch " + str(epoch + 1)
    err = []

    for batch in range(num_batches):
        v_true = dat[:, batch*batch_size:(batch + 1)*batch_size]
        v = v_true

        # apply momentum
        wu_vh *= momentum
        wu_v *= momentum
        wu_h *= momentum

        # positive phase
        h = 1. / (1 + np.exp(-(np.dot(w_vh.T, v) + w_h)))

        wu_vh += np.dot(v, h.T)
        wu_v += v.sum(1)[:, np.newaxis]
        wu_h += h.sum(1)[:, np.newaxis]

        # sample hiddens
        h = 1. * (h > np.random.rand(num_hid, batch_size))

        # negative phase
        v = 1. / (1 + np.exp(-(np.dot(w_vh, h) + w_v)))
        h = 1. / (1 + np.exp(-(np.dot(w_vh.T, v) + w_h)))

        wu_vh -= np.dot(v, h.T)
        wu_v -= v.sum(1)[:, np.newaxis]
        wu_h -= h.sum(1)[:, np.newaxis]

        # update weights
        w_vh += epsilon/batch_size * wu_vh
        w_v += epsilon/batch_size * wu_v
        w_h += epsilon/batch_size * wu_h

        err.append(np.mean((v - v_true)**2))

    print "Mean squared error: " + str(np.mean(err))
    print "Time: " + str(time.time() - start_time)