# HG changeset patch # User Frederic Bastien # Date 1298669913 18000 # Node ID c584d8f8f280c7131e19031ab6277d7b399fad4d # Parent 8f15ef65659810c4c3a614cdd2ab105e3d5bdb6c fixed indentation. diff -r 8f15ef656598 -r c584d8f8f280 pylearn/algorithms/pca_online_estimator.py --- a/pylearn/algorithms/pca_online_estimator.py Fri Feb 25 16:37:48 2011 -0500 +++ b/pylearn/algorithms/pca_online_estimator.py Fri Feb 25 16:38:33 2011 -0500 @@ -7,170 +7,169 @@ # - reevaluation when not at the end of a minibatch class PcaOnlineEstimator(object): - """Online estimation of the leading eigen values/vectors of the covariance - of some samples. + """Online estimation of the leading eigen values/vectors of the covariance + of some samples. - Maintains a moving (with discount) low rank (n_eigen) estimate of the - covariance matrix of some observations. New observations are accumulated - until the batch is complete, at which point the low rank estimate is - reevaluated. + Maintains a moving (with discount) low rank (n_eigen) estimate of the + covariance matrix of some observations. New observations are accumulated + until the batch is complete, at which point the low rank estimate is + reevaluated. - Example: + Example: - pca_esti = pca_online_estimator.PcaOnlineEstimator(dimension_of_the_samples) + pca_esti = pca_online_estimator.PcaOnlineEstimator(dimension_of_the_samples) - for i in range(number_of_samples): - pca_esti.observe(samples[i]) + for i in range(number_of_samples): + pca_esti.observe(samples[i]) - [eigvals, eigvecs] = pca_esti.getLeadingEigen() + [eigvals, eigvecs] = pca_esti.getLeadingEigen() - """ + """ - def __init__(self, n_dim, n_eigen = 10, minibatch_size = 25, gamma = 0.999, regularizer = 1e-6, centering = True): - # dimension of the observations - self.n_dim = n_dim - # rank of the low-rank estimate - self.n_eigen = n_eigen - # how many observations between reevaluations of the low rank estimate - self.minibatch_size = minibatch_size - # the discount factor in the moving estimate - self.gamma = gamma - # regularizer of the covariance estimate - self.regularizer = regularizer - # wether we center the observations or not: obtain leading eigen of - # covariance (centering = True) vs second moment (centering = False) - self.centering = centering + def __init__(self, n_dim, n_eigen = 10, minibatch_size = 25, gamma = 0.999, regularizer = 1e-6, centering = True): + # dimension of the observations + self.n_dim = n_dim + # rank of the low-rank estimate + self.n_eigen = n_eigen + # how many observations between reevaluations of the low rank estimate + self.minibatch_size = minibatch_size + # the discount factor in the moving estimate + self.gamma = gamma + # regularizer of the covariance estimate + self.regularizer = regularizer + # wether we center the observations or not: obtain leading eigen of + # covariance (centering = True) vs second moment (centering = False) + self.centering = centering - # Total number of observations: to compute the normalizer for the mean and - # the covariance. - self.n_observations = 0 - # Index in the current minibatch - self.minibatch_index = 0 + # Total number of observations: to compute the normalizer for the mean and + # the covariance. + self.n_observations = 0 + # Index in the current minibatch + self.minibatch_index = 0 - # Matrix containing on its *rows*: - # - the current unnormalized eigen vector estimates - # - the observations since the last reevaluation - self.Xt = numpy.zeros([self.n_eigen + self.minibatch_size, self.n_dim]) + # Matrix containing on its *rows*: + # - the current unnormalized eigen vector estimates + # - the observations since the last reevaluation + self.Xt = numpy.zeros([self.n_eigen + self.minibatch_size, self.n_dim]) - # The discounted sum of the observations. - self.x_sum = numpy.zeros([self.n_dim]) + # The discounted sum of the observations. + self.x_sum = numpy.zeros([self.n_dim]) - # The Gram matrix of the observations, ie Xt Xt' (since Xt is rowwise) - self.G = numpy.zeros([self.n_eigen + self.minibatch_size, self.n_eigen + self.minibatch_size]) - for i in range(self.n_eigen): - self.G[i,i] = self.regularizer + # The Gram matrix of the observations, ie Xt Xt' (since Xt is rowwise) + self.G = numpy.zeros([self.n_eigen + self.minibatch_size, self.n_eigen + self.minibatch_size]) + for i in range(self.n_eigen): + self.G[i,i] = self.regularizer - # I don't think it's worth "allocating" these 3 next (though they need to be - # declared). I don't know how to do in place operations... + # I don't think it's worth "allocating" these 3 next (though they need to be + # declared). I don't know how to do in place operations... - # Hold the results of the eigendecomposition of the Gram matrix G - # (eigen vectors on columns of V). - self.d = numpy.zeros([self.n_eigen + self.minibatch_size]) - self.V = numpy.zeros([self.n_eigen + self.minibatch_size, self.n_eigen + self.minibatch_size]) + # Hold the results of the eigendecomposition of the Gram matrix G + # (eigen vectors on columns of V). + self.d = numpy.zeros([self.n_eigen + self.minibatch_size]) + self.V = numpy.zeros([self.n_eigen + self.minibatch_size, self.n_eigen + self.minibatch_size]) - # Holds the unnormalized eigenvectors of the covariance matrix before - # they're copied back to Xt. - self.Ut = numpy.zeros([self.n_eigen, self.n_dim]) + # Holds the unnormalized eigenvectors of the covariance matrix before + # they're copied back to Xt. + self.Ut = numpy.zeros([self.n_eigen, self.n_dim]) - def observe(self, x): - assert(numpy.size(x) == self.n_dim) + def observe(self, x): + assert(numpy.size(x) == self.n_dim) - self.n_observations += 1 + self.n_observations += 1 - # Add the *non-centered* observation to Xt. - row = self.n_eigen + self.minibatch_index - self.Xt[row] = x + # Add the *non-centered* observation to Xt. + row = self.n_eigen + self.minibatch_index + self.Xt[row] = x - # Update the discounted sum of the observations. - self.x_sum *= self.gamma - self.x_sum += x + # Update the discounted sum of the observations. + self.x_sum *= self.gamma + self.x_sum += x - # To get the mean, we must normalize the sum by: - # /gamma^(n_observations-1) + /gamma^(n_observations-2) + ... + 1 - normalizer = (1.0 - pow(self.gamma, self.n_observations)) /(1.0 - self.gamma); - #print "normalizer: ", normalizer + # To get the mean, we must normalize the sum by: + # /gamma^(n_observations-1) + /gamma^(n_observations-2) + ... + 1 + normalizer = (1.0 - pow(self.gamma, self.n_observations)) /(1.0 - self.gamma); + #print "normalizer: ", normalizer - # Now center the observation. - # We will lose the first observation as it is the only one in the mean. - if self.centering: - self.Xt[row] -= self.x_sum / normalizer + # Now center the observation. + # We will lose the first observation as it is the only one in the mean. + if self.centering: + self.Xt[row] -= self.x_sum / normalizer - # Multiply the observation by the discount compensator. Basically - # we make this observation look "younger" than the previous ones. The actual - # discount is applied in the reevaluation (and when solving the equations in - # the case of TONGA) by multiplying every direction with the same aging factor. - rn = pow(self.gamma, -0.5*(self.minibatch_index+1)); - self.Xt[row] *= rn + # Multiply the observation by the discount compensator. Basically + # we make this observation look "younger" than the previous ones. The actual + # discount is applied in the reevaluation (and when solving the equations in + # the case of TONGA) by multiplying every direction with the same aging factor. + rn = pow(self.gamma, -0.5*(self.minibatch_index+1)); + self.Xt[row] *= rn - # Update the Gram Matrix. - # The column. - self.G[:row+1,row] = numpy.dot( self.Xt[:row+1,:], self.Xt[row,:].transpose() ) - # The symetric row. - # There are row+1 values, but the diag doesn't need to get copied. - self.G[row,:row] = self.G[:row,row].transpose() + # Update the Gram Matrix. + # The column. + self.G[:row+1,row] = numpy.dot( self.Xt[:row+1,:], self.Xt[row,:].transpose() ) + # The symetric row. + # There are row+1 values, but the diag doesn't need to get copied. + self.G[row,:row] = self.G[:row,row].transpose() - self.minibatch_index += 1 + self.minibatch_index += 1 - if self.minibatch_index == self.minibatch_size: - self.reevaluate() + if self.minibatch_index == self.minibatch_size: + self.reevaluate() - def reevaluate(self): - # TODO do the modifications to handle when this is not true. - assert(self.minibatch_index == self.minibatch_size); + def reevaluate(self): + # TODO do the modifications to handle when this is not true. + assert(self.minibatch_index == self.minibatch_size); - # Regularize - not necessary but in case - for i in range(self.n_eigen + self.minibatch_size): - self.G[i,i] += self.regularizer + # Regularize - not necessary but in case + for i in range(self.n_eigen + self.minibatch_size): + self.G[i,i] += self.regularizer - # The Gram matrix is up to date. Get its low rank eigendecomposition. - # NOTE: the eigenvalues are in ASCENDING order and the vectors are on - # the columns. - # With scipy 0.7, you can ask for only some eigenvalues (the n_eigen top - # ones) but it doesn't look loke it for scipy 0.6. - self.d, self.V = linalg.eigh(self.G) #, overwrite_a=True) + # The Gram matrix is up to date. Get its low rank eigendecomposition. + # NOTE: the eigenvalues are in ASCENDING order and the vectors are on + # the columns. + # With scipy 0.7, you can ask for only some eigenvalues (the n_eigen top + # ones) but it doesn't look loke it for scipy 0.6. + self.d, self.V = linalg.eigh(self.G) #, overwrite_a=True) - # Convert the n_eigen LAST eigenvectors of the Gram matrix contained in V - # into *unnormalized* eigenvectors U of the covariance (unnormalized wrt - # the eigen values, not the moving average). - self.Ut = numpy.dot(self.V[:,-self.n_eigen:].transpose(), self.Xt) + # Convert the n_eigen LAST eigenvectors of the Gram matrix contained in V + # into *unnormalized* eigenvectors U of the covariance (unnormalized wrt + # the eigen values, not the moving average). + self.Ut = numpy.dot(self.V[:,-self.n_eigen:].transpose(), self.Xt) - # Take into account the discount factor. - # Here, minibatch index is minibatch_size. We age everyone. Because of the - # previous multiplications to make some observations "younger" we multiply - # everyone by the same factor. - # TODO VERIFY THIS! - rn = pow(self.gamma, -0.5*(self.minibatch_index+1)) - inv_rn2 = 1.0/(rn*rn) - self.Ut *= 1.0/rn - self.d *= inv_rn2; + # Take into account the discount factor. + # Here, minibatch index is minibatch_size. We age everyone. Because of the + # previous multiplications to make some observations "younger" we multiply + # everyone by the same factor. + # TODO VERIFY THIS! + rn = pow(self.gamma, -0.5*(self.minibatch_index+1)) + inv_rn2 = 1.0/(rn*rn) + self.Ut *= 1.0/rn + self.d *= inv_rn2; - #print "*** Reevaluate! ***" - #normalizer = (1.0 - pow(self.gamma, self.n_observations)) /(1.0 - self.gamma) - #print "normalizer: ", normalizer - #print self.d / normalizer - #print self.Ut # unnormalized eigen vectors (wrt eigenvalues AND moving average). - - # Update Xt, G and minibatch_index - self.Xt[:self.n_eigen,:] = self.Ut + #print "*** Reevaluate! ***" + #normalizer = (1.0 - pow(self.gamma, self.n_observations)) /(1.0 - self.gamma) + #print "normalizer: ", normalizer + #print self.d / normalizer + #print self.Ut # unnormalized eigen vectors (wrt eigenvalues AND moving average). - for i in range(self.n_eigen): - self.G[i,i] = self.d[-self.n_eigen+i] + # Update Xt, G and minibatch_index + self.Xt[:self.n_eigen,:] = self.Ut - self.minibatch_index = 0 + for i in range(self.n_eigen): + self.G[i,i] = self.d[-self.n_eigen+i] + + self.minibatch_index = 0 - # Returns a copy of the current estimate of the eigen values and vectors - # (normalized vectors on rows), normalized by the discounted number of observations. - def getLeadingEigen(self): - # We subtract self.minibatch_index in case this call is not right after a reevaluate call. - normalizer = (1.0 - pow(self.gamma, self.n_observations - self.minibatch_index)) /(1.0 - self.gamma) + # Returns a copy of the current estimate of the eigen values and vectors + # (normalized vectors on rows), normalized by the discounted number of observations. + def getLeadingEigen(self): + # We subtract self.minibatch_index in case this call is not right after a reevaluate call. + normalizer = (1.0 - pow(self.gamma, self.n_observations - self.minibatch_index)) /(1.0 - self.gamma) - eigvals = self.d[-self.n_eigen:] / normalizer - eigvecs = numpy.zeros([self.n_eigen, self.n_dim]) - for i in range(self.n_eigen): - eigvecs[i] = self.Ut[-self.n_eigen+i] / numpy.sqrt(numpy.dot(self.Ut[-self.n_eigen+i], self.Ut[-self.n_eigen+i])) + eigvals = self.d[-self.n_eigen:] / normalizer + eigvecs = numpy.zeros([self.n_eigen, self.n_dim]) + for i in range(self.n_eigen): + eigvecs[i] = self.Ut[-self.n_eigen+i] / numpy.sqrt(numpy.dot(self.Ut[-self.n_eigen+i], self.Ut[-self.n_eigen+i])) - return [eigvals, eigvecs] - + return [eigvals, eigvecs]