Mercurial > ift6266
annotate code_tutoriel/mlp.py @ 34:e61a46996db6
first commit, the module has been tested with testmod.py, it has no timing
| author | ychherawala@smets-w04 |
|---|---|
| date | Tue, 02 Feb 2010 10:00:03 -0500 |
| parents | 827de2cc34f8 |
| children | 4bc5eeec6394 |
| rev | line source |
|---|---|
| 0 | 1 """ |
| 2 This tutorial introduces the multilayer perceptron using Theano. | |
| 3 | |
| 4 A multilayer perceptron is a logistic regressor where | |
| 5 instead of feeding the input to the logistic regression you insert a | |
| 6 intermidiate layer, called the hidden layer, that has a nonlinear | |
| 7 activation function (usually tanh or sigmoid) . One can use many such | |
| 8 hidden layers making the architecture deep. The tutorial will also tackle | |
| 9 the problem of MNIST digit classification. | |
| 10 | |
| 11 .. math:: | |
| 12 | |
| 13 f(x) = G( b^{(2)} + W^{(2)}( s( b^{(1)} + W^{(1)} x))), | |
| 14 | |
| 15 References: | |
| 16 | |
| 17 - textbooks: "Pattern Recognition and Machine Learning" - | |
| 18 Christopher M. Bishop, section 5 | |
| 19 | |
| 20 TODO: recommended preprocessing, lr ranges, regularization ranges (explain | |
| 21 to do lr first, then add regularization) | |
| 22 | |
| 23 """ | |
| 24 __docformat__ = 'restructedtext en' | |
| 25 | |
| 26 | |
| 27 import numpy, cPickle, gzip | |
| 28 | |
| 29 | |
| 30 import theano | |
| 31 import theano.tensor as T | |
| 32 | |
| 33 import time | |
| 34 | |
| 35 import theano.tensor.nnet | |
| 36 | |
| 37 class MLP(object): | |
| 38 """Multi-Layer Perceptron Class | |
| 39 | |
| 40 A multilayer perceptron is a feedforward artificial neural network model | |
| 41 that has one layer or more of hidden units and nonlinear activations. | |
| 42 Intermidiate layers usually have as activation function thanh or the | |
| 43 sigmoid function while the top layer is a softamx layer. | |
| 44 """ | |
| 45 | |
| 46 | |
| 47 | |
| 48 def __init__(self, input, n_in, n_hidden, n_out): | |
| 49 """Initialize the parameters for the multilayer perceptron | |
| 50 | |
| 51 :param input: symbolic variable that describes the input of the | |
| 52 architecture (one minibatch) | |
| 53 | |
| 54 :param n_in: number of input units, the dimension of the space in | |
| 55 which the datapoints lie | |
| 56 | |
| 57 :param n_hidden: number of hidden units | |
| 58 | |
| 59 :param n_out: number of output units, the dimension of the space in | |
| 60 which the labels lie | |
| 61 | |
| 62 """ | |
| 63 | |
| 64 # initialize the parameters theta = (W1,b1,W2,b2) ; note that this | |
| 65 # example contains only one hidden layer, but one can have as many | |
| 66 # layers as he/she wishes, making the network deeper. The only | |
| 67 # problem making the network deep this way is during learning, | |
| 68 # backpropagation being unable to move the network from the starting | |
| 69 # point towards; this is where pre-training helps, giving a good | |
| 70 # starting point for backpropagation, but more about this in the | |
| 71 # other tutorials | |
| 72 | |
| 73 # `W1` is initialized with `W1_values` which is uniformely sampled | |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
74 # from -6./sqrt(n_in+n_hidden) and 6./sqrt(n_in+n_hidden) |
| 0 | 75 # the output of uniform if converted using asarray to dtype |
| 76 # theano.config.floatX so that the code is runable on GPU | |
| 77 W1_values = numpy.asarray( numpy.random.uniform( \ | |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
78 low = -numpy.sqrt(6./(n_in+n_hidden)), \ |
|
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
79 high = numpy.sqrt(6./(n_in+n_hidden)), \ |
| 0 | 80 size = (n_in, n_hidden)), dtype = theano.config.floatX) |
| 81 # `W2` is initialized with `W2_values` which is uniformely sampled | |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
82 # from -6./sqrt(n_hidden+n_out) and 6./sqrt(n_hidden+n_out) |
| 0 | 83 # the output of uniform if converted using asarray to dtype |
| 84 # theano.config.floatX so that the code is runable on GPU | |
| 85 W2_values = numpy.asarray( numpy.random.uniform( | |
|
18
827de2cc34f8
Error on the sign of the lower bound of the initialization of W2
Owner <salahmeister@gmail.com>
parents:
2
diff
changeset
|
86 low = -numpy.sqrt(6./(n_hidden+n_out)), \ |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
87 high= numpy.sqrt(6./(n_hidden+n_out)),\ |
| 0 | 88 size= (n_hidden, n_out)), dtype = theano.config.floatX) |
| 89 | |
| 90 self.W1 = theano.shared( value = W1_values ) | |
| 91 self.b1 = theano.shared( value = numpy.zeros((n_hidden,), | |
| 92 dtype= theano.config.floatX)) | |
| 93 self.W2 = theano.shared( value = W2_values ) | |
| 94 self.b2 = theano.shared( value = numpy.zeros((n_out,), | |
| 95 dtype= theano.config.floatX)) | |
| 96 | |
| 97 # symbolic expression computing the values of the hidden layer | |
| 98 self.hidden = T.tanh(T.dot(input, self.W1)+ self.b1) | |
| 99 | |
| 100 # symbolic expression computing the values of the top layer | |
| 101 self.p_y_given_x= T.nnet.softmax(T.dot(self.hidden, self.W2)+self.b2) | |
| 102 | |
| 103 # compute prediction as class whose probability is maximal in | |
| 104 # symbolic form | |
| 105 self.y_pred = T.argmax( self.p_y_given_x, axis =1) | |
| 106 | |
| 107 # L1 norm ; one regularization option is to enforce L1 norm to | |
| 108 # be small | |
| 109 self.L1 = abs(self.W1).sum() + abs(self.W2).sum() | |
| 110 | |
| 111 # square of L2 norm ; one regularization option is to enforce | |
| 112 # square of L2 norm to be small | |
| 113 self.L2_sqr = (self.W1**2).sum() + (self.W2**2).sum() | |
| 114 | |
| 115 | |
| 116 | |
| 117 def negative_log_likelihood(self, y): | |
| 118 """Return the mean of the negative log-likelihood of the prediction | |
| 119 of this model under a given target distribution. | |
| 120 | |
| 121 .. math:: | |
| 122 | |
| 123 \frac{1}{|\mathcal{D}|}\mathcal{L} (\theta=\{W,b\}, \mathcal{D}) = | |
| 124 \frac{1}{|\mathcal{D}|}\sum_{i=0}^{|\mathcal{D}|} \log(P(Y=y^{(i)}|x^{(i)}, W,b)) \\ | |
| 125 \ell (\theta=\{W,b\}, \mathcal{D}) | |
| 126 | |
| 127 | |
| 128 :param y: corresponds to a vector that gives for each example the | |
| 129 :correct label | |
| 130 """ | |
| 131 return -T.mean(T.log(self.p_y_given_x)[T.arange(y.shape[0]),y]) | |
| 132 | |
| 133 | |
| 134 | |
| 135 | |
| 136 def errors(self, y): | |
| 137 """Return a float representing the number of errors in the minibatch | |
| 138 over the total number of examples of the minibatch | |
| 139 """ | |
| 140 | |
| 141 # check if y has same dimension of y_pred | |
| 142 if y.ndim != self.y_pred.ndim: | |
| 143 raise TypeError('y should have the same shape as self.y_pred', | |
| 144 ('y', target.type, 'y_pred', self.y_pred.type)) | |
| 145 # check if y is of the correct datatype | |
| 146 if y.dtype.startswith('int'): | |
| 147 # the T.neq operator returns a vector of 0s and 1s, where 1 | |
| 148 # represents a mistake in prediction | |
| 149 return T.mean(T.neq(self.y_pred, y)) | |
| 150 else: | |
| 151 raise NotImplementedError() | |
| 152 | |
| 153 | |
| 154 | |
| 155 def sgd_optimization_mnist( learning_rate=0.01, L1_reg = 0.00, \ | |
| 156 L2_reg = 0.0001, n_iter=100): | |
| 157 """ | |
| 158 Demonstrate stochastic gradient descent optimization for a multilayer | |
| 159 perceptron | |
| 160 | |
| 161 This is demonstrated on MNIST. | |
| 162 | |
| 163 :param learning_rate: learning rate used (factor for the stochastic | |
| 164 gradient | |
| 165 | |
| 166 :param L1_reg: L1-norm's weight when added to the cost (see | |
| 167 regularization) | |
| 168 | |
| 169 :param L2_reg: L2-norm's weight when added to the cost (see | |
| 170 regularization) | |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
171 |
|
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
172 :param n_iter: maximal number of iterations ot run the optimizer |
|
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
173 |
|
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
174 """ |
| 0 | 175 |
| 176 # Load the dataset | |
| 177 f = gzip.open('mnist.pkl.gz','rb') | |
| 178 train_set, valid_set, test_set = cPickle.load(f) | |
| 179 f.close() | |
| 180 | |
| 181 # make minibatches of size 20 | |
| 182 batch_size = 20 # sized of the minibatch | |
| 183 | |
| 184 # Dealing with the training set | |
| 185 # get the list of training images (x) and their labels (y) | |
| 186 (train_set_x, train_set_y) = train_set | |
| 187 # initialize the list of training minibatches with empty list | |
| 188 train_batches = [] | |
| 189 for i in xrange(0, len(train_set_x), batch_size): | |
| 190 # add to the list of minibatches the minibatch starting at | |
| 191 # position i, ending at position i+batch_size | |
| 192 # a minibatch is a pair ; the first element of the pair is a list | |
| 193 # of datapoints, the second element is the list of corresponding | |
| 194 # labels | |
| 195 train_batches = train_batches + \ | |
| 196 [(train_set_x[i:i+batch_size], train_set_y[i:i+batch_size])] | |
| 197 | |
| 198 # Dealing with the validation set | |
| 199 (valid_set_x, valid_set_y) = valid_set | |
| 200 # initialize the list of validation minibatches | |
| 201 valid_batches = [] | |
| 202 for i in xrange(0, len(valid_set_x), batch_size): | |
| 203 valid_batches = valid_batches + \ | |
| 204 [(valid_set_x[i:i+batch_size], valid_set_y[i:i+batch_size])] | |
| 205 | |
| 206 # Dealing with the testing set | |
| 207 (test_set_x, test_set_y) = test_set | |
| 208 # initialize the list of testing minibatches | |
| 209 test_batches = [] | |
| 210 for i in xrange(0, len(test_set_x), batch_size): | |
| 211 test_batches = test_batches + \ | |
| 212 [(test_set_x[i:i+batch_size], test_set_y[i:i+batch_size])] | |
| 213 | |
| 214 | |
| 215 ishape = (28,28) # this is the size of MNIST images | |
| 216 | |
| 217 # allocate symbolic variables for the data | |
| 218 x = T.fmatrix() # the data is presented as rasterized images | |
| 219 y = T.lvector() # the labels are presented as 1D vector of | |
| 220 # [long int] labels | |
| 221 | |
| 222 # construct the logistic regression class | |
| 223 classifier = MLP( input=x.reshape((batch_size,28*28)),\ | |
| 224 n_in=28*28, n_hidden = 500, n_out=10) | |
| 225 | |
| 226 # the cost we minimize during training is the negative log likelihood of | |
| 227 # the model plus the regularization terms (L1 and L2); cost is expressed | |
| 228 # here symbolically | |
| 229 cost = classifier.negative_log_likelihood(y) \ | |
| 230 + L1_reg * classifier.L1 \ | |
| 231 + L2_reg * classifier.L2_sqr | |
| 232 | |
| 233 # compiling a theano function that computes the mistakes that are made by | |
| 234 # the model on a minibatch | |
| 235 test_model = theano.function([x,y], classifier.errors(y)) | |
| 236 | |
| 237 # compute the gradient of cost with respect to theta = (W1, b1, W2, b2) | |
| 238 g_W1 = T.grad(cost, classifier.W1) | |
| 239 g_b1 = T.grad(cost, classifier.b1) | |
| 240 g_W2 = T.grad(cost, classifier.W2) | |
| 241 g_b2 = T.grad(cost, classifier.b2) | |
| 242 | |
| 243 # specify how to update the parameters of the model as a dictionary | |
| 244 updates = \ | |
| 245 { classifier.W1: classifier.W1 - learning_rate*g_W1 \ | |
| 246 , classifier.b1: classifier.b1 - learning_rate*g_b1 \ | |
| 247 , classifier.W2: classifier.W2 - learning_rate*g_W2 \ | |
| 248 , classifier.b2: classifier.b2 - learning_rate*g_b2 } | |
| 249 | |
| 250 # compiling a theano function `train_model` that returns the cost, but in | |
| 251 # the same time updates the parameter of the model based on the rules | |
| 252 # defined in `updates` | |
| 253 train_model = theano.function([x, y], cost, updates = updates ) | |
| 254 n_minibatches = len(train_batches) | |
| 255 | |
| 256 # early-stopping parameters | |
| 257 patience = 10000 # look as this many examples regardless | |
| 258 patience_increase = 2 # wait this much longer when a new best is | |
| 259 # found | |
| 260 improvement_threshold = 0.995 # a relative improvement of this much is | |
| 261 # considered significant | |
| 262 validation_frequency = n_minibatches # go through this many | |
| 263 # minibatche before checking the network | |
| 264 # on the validation set; in this case we | |
| 265 # check every epoch | |
| 266 | |
| 267 | |
| 268 best_params = None | |
| 269 best_validation_loss = float('inf') | |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
270 best_iter = 0 |
| 0 | 271 test_score = 0. |
| 272 start_time = time.clock() | |
| 273 # have a maximum of `n_iter` iterations through the entire dataset | |
| 274 for iter in xrange(n_iter* n_minibatches): | |
| 275 | |
| 276 # get epoch and minibatch index | |
| 277 epoch = iter / n_minibatches | |
| 278 minibatch_index = iter % n_minibatches | |
| 279 | |
| 280 # get the minibatches corresponding to `iter` modulo | |
| 281 # `len(train_batches)` | |
| 282 x,y = train_batches[ minibatch_index ] | |
| 283 cost_ij = train_model(x,y) | |
| 284 | |
| 285 if (iter+1) % validation_frequency == 0: | |
| 286 # compute zero-one loss on validation set | |
| 287 this_validation_loss = 0. | |
| 288 for x,y in valid_batches: | |
| 289 # sum up the errors for each minibatch | |
| 290 this_validation_loss += test_model(x,y) | |
| 291 # get the average by dividing with the number of minibatches | |
| 292 this_validation_loss /= len(valid_batches) | |
| 293 | |
| 294 print('epoch %i, minibatch %i/%i, validation error %f %%' % \ | |
| 295 (epoch, minibatch_index+1, n_minibatches, \ | |
| 296 this_validation_loss*100.)) | |
| 297 | |
| 298 | |
| 299 # if we got the best validation score until now | |
| 300 if this_validation_loss < best_validation_loss: | |
| 301 | |
| 302 #improve patience if loss improvement is good enough | |
| 303 if this_validation_loss < best_validation_loss * \ | |
| 304 improvement_threshold : | |
| 305 patience = max(patience, iter * patience_increase) | |
| 306 | |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
307 # save best validation score and iteration number |
| 0 | 308 best_validation_loss = this_validation_loss |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
309 best_iter = iter |
|
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
310 |
| 0 | 311 # test it on the test set |
| 312 test_score = 0. | |
| 313 for x,y in test_batches: | |
| 314 test_score += test_model(x,y) | |
| 315 test_score /= len(test_batches) | |
| 316 print((' epoch %i, minibatch %i/%i, test error of best ' | |
| 317 'model %f %%') % | |
| 318 (epoch, minibatch_index+1, n_minibatches, | |
| 319 test_score*100.)) | |
| 320 | |
| 321 if patience <= iter : | |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
322 break |
| 0 | 323 |
| 324 end_time = time.clock() | |
|
2
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
325 print(('Optimization complete. Best validation score of %f %% ' |
|
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
326 'obtained at iteration %i, with test performance %f %%') % |
|
bcc87d3e33a3
adding latest tutorial code
Dumitru Erhan <dumitru.erhan@gmail.com>
parents:
0
diff
changeset
|
327 (best_validation_loss * 100., best_iter, test_score*100.)) |
| 0 | 328 print ('The code ran for %f minutes' % ((end_time-start_time)/60.)) |
| 329 | |
| 330 | |
| 331 if __name__ == '__main__': | |
| 332 sgd_optimization_mnist() | |
| 333 |