Mercurial > pylearn

comparison kernel_regression.py @ 421:e01f17be270a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Kernel regression learning algorithm
author Yoshua Bengio <bengioy@iro.umontreal.ca>
date Sat, 19 Jul 2008 10:11:22 -0400
parents
children 32c5f87bc54e
comparison
equal deleted inserted replaced
420:040cb796f4e0 421:e01f17be270a
1 """
2 Implementation of kernel regression:
3 """
4
5 from pylearn.learner import OfflineLearningAlgorithm
6 from theano import tensor as T
7 from nnet_ops import prepend_1_to_each_row
8 from theano.scalar import as_scalar
9 from common.autoname import AutoName
10 import theano
11 import numpy
12
13 class KernelRegression(OfflineLearningAlgorithm):
14 """
15 Implementation of kernel regression:
16 * the data are n (x_t,y_t) pairs and we want to estimate E[y|x]
17 * the predictor computes
18 f(x) = b + \sum_{t=1}^n \alpha_t K(x,x_t)
19 with free parameters b and alpha, training inputs x_t,
20 and kernel function K (gaussian by default).
21 Clearly, each prediction involves O(n) computations.
22 * the learner chooses b and alpha to minimize
23 lambda alpha' G' G alpha + \sum_{t=1}^n (f(x_t)-y_t)^2
24 where G is the matrix with entries G_ij = K(x_i,x_j).
25 The first (L2 regularization) term is the squared L2
26 norm of the primal weights w = \sum_t \alpha_t phi(x_t)
27 where phi is the function s.t. K(u,v)=phi(u).phi(v).
28 * this involves solving a linear system with (n+1,n+1)
29 matrix, which is an O(n^3) computation. In addition,
30 that linear system matrix requires O(n^2) memory.
31 So this learning algorithm should be used only for
32 small datasets.
33 * the linear system is
34 (M + lambda I_n) theta = (1, y)'
35 where theta = (b, alpha), I_n is the (n+1)x(n+1) matrix that is the identity
36 except with a 0 at (0,0), M is the matrix with G in the sub-matrix starting
37 at (1,1), 1's in column 0, except for a value of n at (0,0), and sum_i G_{i,j}
38 in the rest of row 0.
39
40 Note that this is gives an estimate of E[y|x,training_set] that is the
41 same as obtained with a Gaussian process regression. The GP
42 regression would also provide a Bayesian Var[y|x,training_set].
43 It corresponds to an assumption that f is a random variable
44 with Gaussian (process) prior distribution with covariance
45 function K. Because we assume Gaussian noise we obtain a Gaussian
46 posterior for f (whose mean is computed here).
47
48
49 Usage:
50
51 kernel_regressor=KernelRegression(L2_regularizer=0.1,kernel=GaussianKernel(gamma=0.5))
52 kernel_predictor=kernel_regressor(training_set)
53 all_results_dataset=kernel_predictor(test_set) # creates a dataset with "output" and "squared_error" field
54 outputs = kernel_predictor.compute_outputs(inputs) # inputs and outputs are numpy arrays
55 outputs, errors = kernel_predictor.compute_outputs_and_errors(inputs,targets)
56 errors = kernel_predictor.compute_errors(inputs,targets)
57 mse = kernel_predictor.compute_mse(inputs,targets)
58
59
60
61 The training_set must have fields "input" and "target".
62 The test_set must have field "input", and needs "target" if
63 we want to compute the squared errors.
64
65 The predictor parameters are obtained analytically from the training set.
66 Training is only done on a whole training set rather than on minibatches
67 (no online implementation).
68
69 The dataset fields expected and produced by the learning algorithm and the trained model
70 are the following:
71
72 - Input and output dataset fields (example-wise quantities):
73
74 - 'input' (always expected as an input_dataset field)
75 - 'target' (always expected by the learning algorithm, optional for learned model)
76 - 'output' (always produced by learned model)
77 - 'squared_error' (optionally produced by learned model if 'target' is provided)
78 = example-wise squared error
79 """
80 def __init__(self, kernel=None, L2_regularizer=0, gamma=1):
81 self.kernel = kernel
82 self.L2_regularizer=L2_regularizer
83 self.gamma = gamma # until we fix things, the kernel type is fixed, Gaussian
84 self.equations = KernelRegressionEquations()
85
86 def __call__(self,trainset):
87 n_examples = len(trainset)
88 first_example = trainset[0]
89 n_inputs = first_example['input'].size
90 n_outputs = first_example['target'].size
91 M = numpy.zeros((n_examples+1,n_examples+1))
92 Y = numpy.zeros((n_examples+1,n_outputs))
93 for i in xrange(n_inputs):
94 M[i+1,i+1]=self.L2_regularizer
95 data = trainset.fields()
96 train_inputs = numpy.array(data['input'])
97 Y[0]=1
98 Y[1:,:] = numpy.array(data['target'])
99 M,train_inputs_square=self.equations.compute_system_matrix(train_inputs,M)
100 theta=numpy.linalg.solve(M,Y)
101 return KernelPredictor(theta,self.gamma, train_inputs, train_inputs_square)
102
103 class KernelPredictorEquations(AutoName):
104 train_inputs = T.matrix() # n_examples x n_inputs
105 train_inputs_square = T.vector() # n_examples
106 inputs = T.matrix() # minibatchsize x n_inputs
107 targets = T.matrix() # minibatchsize x n_outputs
108 theta = T.matrix() # (n_examples+1) x n_outputs
109 gamma = T.scalar()
110 inv_gamma2 = 1./(gamma*gamma)
111 b = theta[0]
112 alpha = theta[1:,:]
113 inputs_square = T.sum(inputs*inputs,axis=1)
114 Kx = exp(-(train_inputs_square-2*dot(inputs,train_inputs.T)+inputs_square)*inv_gamma2)
115 outputs = T.dot(Kx,alpha) + b # minibatchsize x n_outputs
116 squared_errors = T.sum(T.sqr(targets-outputs),axis=1)
117
118 __compiled = False
119 @classmethod
120 def compile(cls,linker='c|py'):
121 if cls.__compiled:
122 return
123 def fn(input_vars,output_vars):
124 return staticmethod(theano.function(input_vars,output_vars, linker=linker))
125
126 cls.compute_outputs = fn([cls.inputs,cls.theta,cls.gamma,cls.train_inputs,cls.train_inputs_square],[cls.outputs])
127 cls.compute_errors = fn([cls.outputs,cls.targets],[cls.squared_errors])
128
129 cls.__compiled = True
130
131 def __init__(self):
132 self.compile()
133
134 class KernelRegressionEquations(KernelPredictorEquations):
135 # P = KernelPredictorEquations
136 M = T.matrix() # (n_examples+1) x (n_examples+1)
137 inputs = T.matrix() # n_examples x n_inputs
138 G = M[1:,1:]
139 new_G = gemm(G,1.,inputs,inputs.T,1.)
140 M2 = T.add_inplace(M,new_G)
141 M2[0,0] = M.shape[0]
142 M2[1:,0] = 1
143 M2[0,1:] = T.sum(G,axis=0)
144 inputs_square = T.sum(inputs*inputs,axis=1)
145
146 __compiled = False
147
148 @classmethod
149 def compile(cls,linker='c|py'):
150 if cls.__compiled:
151 return
152 def fn(input_vars,output_vars):
153 return staticmethod(theano.function(input_vars,output_vars, linker=linker))
154
155 cls.compute_system_matrix = fn([cls.inputs,cls.M],[cls.M2,cls.inputs_square])
156
157 cls.__compiled = True
158
159 def __init__(self):
160 self.compile()
161
162 class KernelPredictor(object):
163 """
164 A kernel predictor has parameters theta (a bias vector and a weight matrix alpha)
165 it can use to make a non-linear prediction (according to the KernelPredictorEquations).
166 It can compute its output (bias + alpha * kernel(train_inputs,input) and a squared error (||output - target||^2).
167 """
168 def __init__(self, theta, gamma, train_inputs, train_inputs_square):
169 self.theta=theta
170 self.gamma=gamma
171 self.train_inputs=train_inputs
172 self.train_inputs_square=train_inputs_square
173 self.equations = LinearPredictorEquations()
174
175 def compute_outputs(self,inputs):
176 return self.equations.compute_outputs(inputs,self.theta,self.gamma,self.train_inputs,self.train_inputs_square)
177 def compute_errors(self,inputs,targets):
178 return self.equations.compute_errors(self.compute_outputs(inputs),targets)
179 def compute_outputs_and_errors(self,inputs,targets):
180 outputs = self.compute_outputs(inputs)
181 return [outputs,self.equations.compute_errors(outputs,targets)]
182 def compute_mse(self,inputs,targets):
183 errors = self.compute_errors(inputs,targets)
184 return numpy.sum(errors)/errors.size
185
186 def __call__(self,dataset,output_fieldnames=None,cached_output_dataset=False):
187 assert dataset.hasFields(["input"])
188 if output_fieldnames is None:
189 if dataset.hasFields(["target"]):
190 output_fieldnames = ["output","squared_error"]
191 else:
192 output_fieldnames = ["output"]
193 output_fieldnames.sort()
194 if output_fieldnames == ["squared_error"]:
195 f = self.compute_errors
196 elif output_fieldnames == ["output"]:
197 f = self.compute_outputs
198 elif output_fieldnames == ["output","squared_error"]:
199 f = self.compute_outputs_and_errors
200 else:
201 raise ValueError("unknown field(s) in output_fieldnames: "+str(output_fieldnames))
202
203 ds=ApplyFunctionDataSet(dataset,f,output_fieldnames)
204 if cached_output_dataset:
205 return CachedDataSet(ds)
206 else:
207 return ds
208
209