Mercurial > ift6266
annotate code_tutoriel/mlp.py @ 160:68160fd149fe
Simple script to doctest all modules beneath ift6266.
This ecxcludes scripts since most of the files launch computations
when imported and some even start a lot of jobs on the cluster (This
happened to me when testing).
author | Arnaud Bergeron <abergeron@gmail.com> |
---|---|
date | Thu, 25 Feb 2010 17:27:49 -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 |