- committed local changes

git-svn-id: http://moon:8086/svn/software/trunk/projects/RBM@270 b431acfa-c32f-4a4a-93f1-934dc6c82436
This commit is contained in:
2015-05-28 17:41:06 +00:00
parent 549089a440
commit 7fada2a227
10 changed files with 467 additions and 640 deletions
+180 -334
View File
@@ -19,6 +19,8 @@ using namespace Eigen;
void mylog(const char* format, ...);
#define printf mylog
#define EPSILON_SIGMA 0.05
class Rbm;
class RbmListener
@@ -53,9 +55,18 @@ public:
, m_useProbsForHiddenReconstruction(false)
, m_doSparse(false)
, m_doNormalizeData(false)
, m_doLearnVariance(false)
, m_numGibbs(1)
{
Noise_Init(&m_noise, 0x32727155);
VectorXd a(4);
a << 1, 2, 3, 4;
VectorXd b(4);
b.array() = -a.array().exp();
cout << b << endl;
}
~Rbm()
@@ -64,173 +75,6 @@ public:
Noise_Free(&m_noise);
}
void train(const LayerArray<VisibleLayer> &batch, uint32_t numEpochs, double sigmaMin = 0.05)
{
uint32_t t, i;
uint32_t epoch;
uint32_t gibbs;
double sigma;
VisibleLayer v(m_w.getNumVisible());
HiddenLayer h(m_w.getNumHidden());
VectorXd sumBiasV(m_w.getNumVisible());
VectorXd deltaBiasV(m_w.getNumVisible());
VectorXd sumBiasH(m_w.getNumHidden());
VectorXd deltaBiasH(m_w.getNumHidden());
MatrixXd sumWeights(m_w.getNumVisible(), m_w.getNumHidden());
MatrixXd deltaWeights(m_w.getNumVisible(), m_w.getNumHidden());
MatrixXd diffErr(1, m_w.getNumVisible());
const LayerArray<VisibleLayer> &vt = batch;
sigma = m_sigma;
double dProgress = 1.0/numEpochs;
double kTrain = 1.0/vt.getSize();
m_progress = 0;
deltaWeights.fill(0);
deltaBiasV.fill(0);
deltaBiasH.fill(0);
m_doCancel = false;
for (epoch=0; epoch < numEpochs; epoch++)
{
double err = 0;
if (m_doCancel)
{
m_doCancel = false;
break;
}
sumWeights.fill(0);
sumBiasV.fill(0);
sumBiasH.fill(0);
for (i=0; i < vt.getSize(); i++)
{
t = i;
h.probsUpdateLogistic(vt[t], m_w, m_lambda, sigma);
// Create hidden layer base on training data
if (m_doRaoBlackwell)
{
h.states() = h.probs();
}
else
{
h.statesUpdateStochastic();
}
// Update weights (positive phase)
sumWeights += vt[t].states() * h.states().transpose();
sumBiasV += vt[t].states();
sumBiasH += h.states();
diffErr = vt[t].states();
for (gibbs=0; gibbs < m_numGibbs; gibbs++)
{
h.statesUpdateStochastic();
// Create visible reconstruction (a fantasy...)
if (m_useProbsForHiddenReconstruction)
{
if (m_useVisibleGaussian)
{
v.probsUpdateGaussian(h, m_w, m_lambda, sigma);
}
else
{
v.probsUpdateLogistic(h, m_w, m_lambda, sigma);
}
v.states() = v.probs();
}
else
{
if (m_useVisibleGaussian)
{
v.sampleGaussian(h, m_w, m_lambda, sigma);
}
else
{
v.probsUpdateLogistic(h, m_w, m_lambda, sigma);
v.statesUpdateStochastic();
}
}
// Create hidden reconstruction
h.probsUpdateLogistic(v, m_w, m_lambda, sigma);
}
// Update weights (negative phase)
if (m_doRaoBlackwell)
{
h.states() = h.probs();
}
else
{
h.statesUpdateStochastic();
}
sumWeights -= v.states() * h.states().transpose();
sumBiasV -= v.states();
sumBiasH -= h.states();
diffErr -= v.states();
diffErr.array() *= diffErr.array();
err += diffErr.sum();
} // TrainingSize
deltaWeights = m_momentum*deltaWeights + m_muWeights*(kTrain*sumWeights - m_weightDecay*m_w.weights());
m_w.weights() += deltaWeights;
deltaBiasV = m_momentum*deltaBiasV + m_muWeights*kTrain*sumBiasV;
m_w.visibleBias() += deltaBiasV;
if (m_doSparse)
{
HiddenLayer th(m_w.getNumHidden());
VectorXd m(m_w.getNumHidden());
m.fill(0);
for (i=0; i < vt.getSize(); i++)
{
th.probsUpdateLogistic(vt[i], m_w, m_lambda, sigma);
m += th.probs();
}
m /= i;
sumBiasH = m_sparsity - m.array();
deltaBiasH = m_momentum*deltaBiasH + m_muSparsity*sumBiasH;
cout << "Mean(" << m_sparsity << ") = " << (double)m.array().mean() << endl;
cout << m << endl;
}
else
{
deltaBiasH = m_momentum*deltaBiasH + m_muWeights*kTrain*sumBiasH;
}
m_w.hiddenBias() += deltaBiasH;
if (sigma > sigmaMin)
{
sigma *= m_sigmaDecay;
}
m_progress += dProgress;
if (m_pListener)
{
m_pListener->onEpochTrained(*this);
}
cout << "err =" << endl;
cout << err << endl;
} // Number of epochs
}
void sample(MatrixXd &src)
{
uint32_t i;
@@ -241,61 +85,144 @@ public:
}
}
void probsLogistic(MatrixXd &src, double lambda, double sigma)
void probsLogistic(MatrixXd &src)
{
double var = sigma*sigma;
src.array() *= -lambda/var;
src.array() = src.array().exp();
src.array() = (-src.array()).exp();
src.array() += 1;
src.array() = 1.0/src.array();
}
void sampleGaussian(MatrixXd &src, double lambda, double sigma)
void probsLogistic(RowVectorXd &src)
{
src.array() = (-src.array()).exp();
src.array() += 1;
src.array() = 1.0/src.array();
}
void probsLogistic(MatrixXd &src, const MatrixXd &sigma)
{
src.array() /= (sigma.array() + EPSILON_SIGMA);
src.array() = (-src.array()).exp();
src.array() += 1;
src.array() = 1.0/src.array();
}
void probsLogistic(RowVectorXd &src, const RowVectorXd &sigma)
{
src.array() /= (sigma.array() + EPSILON_SIGMA);
src.array() = (-src.array()).exp();
src.array() += 1;
src.array() = 1.0/src.array();
}
void probsGaussian(MatrixXd &src, const MatrixXd &sigma)
{
src.array() = 1 - src.array();
src.array() *= src.array();
src.array() *= -0.5;
MatrixXd var = sigma;
var.array() += EPSILON_SIGMA;
var.array() *= var.array();
src.array() /= var.array();
src.array() = src.array().exp();
MatrixXd k = var;
k.array() *= 2*3.14159265359;
k.array() = k.array().sqrt();
k.array() = 1.0/k.array();
src.array() *= k.array();
}
void probsGaussian(RowVectorXd &src, const RowVectorXd &sigma)
{
src.array() = 1 - src.array();
src.array() *= src.array();
src.array() *= -0.5;
RowVectorXd var = sigma;
var.array() += EPSILON_SIGMA;
var.array() *= var.array();
src.array() /= var.array();
src.array() = src.array().exp();
RowVectorXd k = var;
k.array() *= 2*3.14159265359;
k.array() = k.array().sqrt();
k.array() = 1.0/k.array();
src.array() *= k.array();
}
void sampleGaussian(MatrixXd &src, const MatrixXd &sigma)
{
uint32_t i;
for (i=0; i < src.array().size(); i++)
{
src.array()(i) = sigma*Noise_Gaussian(&m_noise) + lambda*src.array()(i);
src.array()(i) = sigma(i)*Noise_Gaussian(&m_noise) + src.array()(i);
}
}
void normalizeData(MatrixXd &src, double mu, double sigma)
RowVectorXd normalizeData(RowVectorXd const &src, RowVectorXd const &mu, RowVectorXd const &var)
{
uint32_t i;
uint32_t size = src.rows();
double mean;
double stdDev;
// Remove mean
RowVectorXd res = src - mu;
// res.array() /= var.array() + EPSILON_SIGMA;
for (i=0; i < size; i++)
{
mean = src.row(i).array().mean();
src.row(i).array() -= mean;
src.row(i).array() += mu;
}
for (i=0; i < size; i++)
{
src.row(i).array() *= src.row(i).array();
}
for (i=0; i < size; i++)
{
stdDev = sqrt(src.row(i).array().mean());
src.row(i).array() /= stdDev;
src.row(i).array() *= sigma;
}
// cout << __PRETTY_FUNCTION__ << ": " << res << endl;
return res;
}
void train2(const LayerArray<VisibleLayer> &vt, uint32_t numEpochs, uint32_t batchSize, double sigmaMin = 0.05)
RowVectorXd calcMean(MatrixXd const &batch)
{
// Remove mean
RowVectorXd res = batch.colwise().mean();
// cout << __PRETTY_FUNCTION__ << ": " << res << endl;
return res;
}
RowVectorXd calcSigma(MatrixXd const &batch)
{
MatrixXd x = batch.rowwise() - batch.colwise().mean();
x.array() *= x.array();
RowVectorXd res = x.colwise().mean().array().sqrt();
// cout << __PRETTY_FUNCTION__ << ": " << res << endl;
return res;
}
MatrixXd calcZ(MatrixXd &v, MatrixXd &h)
{
MatrixXd t1(v.rows(), m_w.getNumVisible());
t1 = v - m_w.visibleBias().transpose().replicate(v.rows(), 1);
t1.array() *= t1.array();
t1.array() *= 0.5;
t1 -= (h * m_w.weights().transpose());
return t1;
}
void train(const LayerArray &vt, uint32_t numEpochs, uint32_t batchSize, double sigmaMin = 0.05)
{
uint32_t t, i;
uint32_t epoch;
uint32_t gibbs;
double sigma = m_sigma;
double dProgress = 1.0/numEpochs;
double kTrain = 1.0/vt.getSize();
// if (batchSize > vt.getSize())
batchSize = vt.getSize();
MatrixXd v(batchSize, m_w.getNumVisible());
@@ -306,6 +233,7 @@ public:
MatrixXd sumBiasH(1, m_w.getNumHidden());
MatrixXd sumWeights(m_w.getNumVisible(), m_w.getNumHidden());
MatrixXd deltaVar(MatrixXd::Zero(1, m_w.getNumVisible()));
MatrixXd deltaBiasV(MatrixXd::Zero(1, m_w.getNumVisible()));
MatrixXd deltaBiasH(MatrixXd::Zero(1, m_w.getNumHidden()));
MatrixXd deltaWeights(MatrixXd::Zero(m_w.getNumVisible(), m_w.getNumHidden()));
@@ -315,15 +243,22 @@ public:
m_progress = 0;
m_doCancel = false;
for (i=0; i < batchSize; i++)
{
// t = (uint32_t)(0.5 + (vt.getSize()-1)*Noise_Uniform(&m_noise));
batch.row(i) = vt[i].states();
batch = vt.data();
m_w.mean() = calcMean(batch);
if (m_doLearnVariance)
{
m_w.sigma() = calcSigma(batch);
}
if (m_doNormalizeData)
{
normalizeData(batch, 0.0, m_sigma);
for (i=0; i < batchSize; i++)
{
RowVectorXd x = batch.row(i);
batch.row(i) = normalizeData(x, m_w.mean(), m_w.sigma());
}
}
for (epoch=0; epoch < numEpochs; epoch++)
@@ -339,13 +274,14 @@ public:
// Create hidden layer base on training data
h = v * m_w.weights();
h += m_w.hiddenBias().transpose().replicate(batchSize, 1);
probsLogistic(h, m_lambda, sigma);
h += m_w.hiddenBias().replicate(batchSize, 1);
probsLogistic(h);
if (!m_doRaoBlackwell)
{
sample(h);
}
// Update weights (positive phase)
sumBiasV = v.colwise().sum();
if (!m_doSparse)
@@ -361,15 +297,18 @@ public:
// Create visible reconstruction (a fantasy...)
v = h * m_w.weights().transpose();
v += m_w.visibleBias().transpose().replicate(batchSize, 1);
v += m_w.visibleBias().replicate(batchSize, 1);
if (m_useVisibleGaussian)
{
sampleGaussian(v, m_lambda, sigma);
if (!m_useProbsForHiddenReconstruction)
{
sampleGaussian(v, m_w.sigma().replicate(batchSize, 1));
}
}
else
{
probsLogistic(v, m_lambda, sigma);
probsLogistic(v, m_w.sigma().replicate(batchSize, 1));
if (!m_useProbsForHiddenReconstruction)
{
sample(v);
@@ -378,8 +317,8 @@ public:
// Create hidden reconstruction
h = v * m_w.weights();
h += m_w.hiddenBias().transpose().replicate(batchSize, 1);
probsLogistic(h, m_lambda, sigma);
h += m_w.hiddenBias().replicate(batchSize, 1);
probsLogistic(h);
}
if (!m_doRaoBlackwell)
@@ -404,8 +343,8 @@ public:
if (m_doSparse)
{
h = v * m_w.weights();
h += m_w.hiddenBias().transpose().replicate(batchSize, 1);
probsLogistic(h, m_lambda, sigma);
h += m_w.hiddenBias().replicate(batchSize, 1);
probsLogistic(h);
sumBiasH.fill(m_sparsity);
sumBiasH -= h.colwise().mean();
@@ -421,9 +360,9 @@ public:
}
m_w.hiddenBias() += deltaBiasH;
if (sigma > sigmaMin)
if (m_w.sigma()[0] > sigmaMin)
{
sigma *= m_sigmaDecay;
m_w.sigma().array() *= m_sigmaDecay;
}
m_progress += dProgress;
@@ -450,148 +389,45 @@ public:
{
double energy;
energy = m_w.visibleBias().transpose() * visible;
energy += m_w.hiddenBias().transpose() * hidden;
energy = m_w.visibleBias() * visible;
energy += m_w.hiddenBias() * hidden;
energy += visible.transpose() * m_w.weights() * hidden;
return -energy/(m_sigma*m_sigma);
}
void prob(LayerArray<VisibleLayer> &vts)
RowVectorXd toHidden(const RowVectorXd& v)
{
uint32_t i, j;
double z;
double p;
RowVectorXd h(m_w.getNumHidden());
RowVectorXd vn(m_w.getNumVisible());
HiddenLayer *h = new HiddenLayer[vts.getSize()];
h = v * m_w.weights();
h += m_w.hiddenBias();
probsLogistic(h);
// Create hidden layer activations based on training data
for (j=0; j < vts.getSize(); j++)
{
h[j].setNumUnits(m_w.getNumHidden());
h[j].probsUpdateLogistic(vts.getAt(j), m_w, m_lambda, m_sigma);
// h[j].statesAssignfromProbs();
h[j].statesUpdateStochastic();
}
printf("pi(t) = (pi^, v>)\n");
for (j=0; j < vts.getSize(); j++)
{
cout << h[j].probs() << endl;
}
cout << endl;
printf("si(t) = (si^, v>)\n");
for (j=0; j < vts.getSize(); j++)
{
cout << h[j].states() << endl;
}
cout << endl;
printf("p(v) = (t^, v>)\n");
for (i=0; i < vts.getSize(); i++)
{
z = 0;
for (j=0; j < vts.getSize(); j++)
{
z += exp(-getEnergy(vts.getAt(j).states(), h[i].states()));
}
for (j=0; j < vts.getSize(); j++)
{
p = exp(-getEnergy(vts.getAt(j).states(), h[i].states()))/z;
cout << p << endl;
}
cout << endl;
}
cout << endl;
// Reconstruct
for (i=0; i < vts.getSize(); i++)
{
vts.getAt(i).probsUpdateLogistic(h[i], m_w, m_lambda, m_sigma);
}
printf("A fantasy... (v^, t>)\n");
for (j=0; j < vts.getSize(); j++)
{
cout << vts.getAt(j).probs() << endl;
}
delete [] h;
return h;
}
VectorXd toHidden(const VectorXd& visible)
RowVectorXd toVisible(const RowVectorXd& h)
{
HiddenLayer th(m_w.getNumHidden());
VisibleLayer tv(m_w.getNumVisible(), (const VectorXd*)&visible);
RowVectorXd v(m_w.getNumVisible());
v = h * m_w.weights().transpose();
v += m_w.visibleBias();
th.probsUpdateLogistic(tv, m_w, m_lambda, m_sigma);
return th.probs();
}
VectorXd toVisible(const VectorXd& hidden)
{
HiddenLayer th(m_w.getNumHidden(), (const VectorXd*)&hidden);
VisibleLayer tv(m_w.getNumVisible());
if (m_useVisibleGaussian)
{
tv.probsUpdateGaussian(th, m_w, m_lambda, m_sigma);
// probsGaussian(v, m_w.sigma());
}
else
{
tv.probsUpdateLogistic(th, m_w, m_lambda, m_sigma);
probsLogistic(v, m_w.sigma());
}
return tv.probs();
}
VectorXd expectHidden(VectorXd visible, uint32_t numIter)
{
uint32_t i;
VisibleLayer v(m_w.getNumVisible(), (const VectorXd*)&visible);
HiddenLayer h(m_w.getNumHidden());
for (i=0; i < numIter; i++)
{
h.probsUpdateLogistic(v, (Weights&)m_w, m_lambda, m_sigma);
if (m_useVisibleGaussian)
{
v.probsUpdateGaussian(h, (Weights&)m_w, m_lambda, m_sigma);
}
else
{
v.probsUpdateLogistic(h, (Weights&)m_w, m_lambda, m_sigma);
}
}
return h.probs();
}
VectorXd expectVisible(VectorXd visible, uint32_t numIter)
{
uint32_t i;
VisibleLayer v(m_w.getNumVisible(), (const VectorXd*)&visible);
HiddenLayer h(m_w.getNumHidden());
for (i=0; i < numIter; i++)
{
h.probsUpdateLogistic(v, (Weights&)m_w, m_lambda, m_sigma);
if (m_useVisibleGaussian)
{
v.probsUpdateGaussian(h, (Weights&)m_w, m_lambda, m_sigma);
}
else
{
v.probsUpdateLogistic(h, (Weights&)m_w, m_lambda, m_sigma);
}
}
return v.probs();
return v;
}
void setSigma(double value)
{
m_sigma = value;
m_sigma = value;
}
void setSigmaDecay(double value)
@@ -639,6 +475,15 @@ public:
m_doNormalizeData = flag;
}
void setDoLearnVariance(bool flag)
{
m_doLearnVariance = flag;
if (!flag)
{
m_w.sigma().fill(m_sigma);
}
}
void setNumGibbs(uint32_t value)
{
m_numGibbs = value;
@@ -683,6 +528,7 @@ private:
bool m_useProbsForHiddenReconstruction;
bool m_doSparse;
bool m_doNormalizeData;
bool m_doLearnVariance;
volatile bool m_doCancel;
uint32_t m_numGibbs;