jens
93fa8c64ce
- Biases are initialized with zero - developing full matrix calculation in train2() git-svn-id: http://moon:8086/svn/software/trunk/projects/RBM@43 b431acfa-c32f-4a4a-93f1-934dc6c82436
611 lines
12 KiB
C++
611 lines
12 KiB
C++
/*
|
|
* Rbm.hpp
|
|
*
|
|
* Created on: 21.09.2014
|
|
* Author: jens
|
|
*/
|
|
|
|
#ifndef RBM_HPP_
|
|
#define RBM_HPP_
|
|
|
|
#include "VisibleLayer.hpp"
|
|
#include "HiddenLayer.hpp"
|
|
#include "Weights.hpp"
|
|
#include <cmath>
|
|
#include <Eigen/Dense>
|
|
|
|
using namespace Eigen;
|
|
|
|
void mylog(const char* format, ...);
|
|
#define printf mylog
|
|
|
|
class Rbm;
|
|
|
|
class RbmListener
|
|
{
|
|
public:
|
|
RbmListener() {}
|
|
virtual ~RbmListener()
|
|
{
|
|
}
|
|
|
|
virtual void onEpochTrained(const Rbm &obj) = 0;
|
|
};
|
|
|
|
class Rbm
|
|
{
|
|
public:
|
|
Rbm(Weights &weights, RbmListener *pListener = nullptr)
|
|
: m_w(weights)
|
|
, m_pListener(pListener)
|
|
, m_progress(0)
|
|
, m_sigma(1.0)
|
|
, m_sigmaDecay(1.0)
|
|
, m_weightDecay(0.0)
|
|
, m_lambda(1.0)
|
|
, m_sparsity(0)
|
|
, m_muWeights(0.01)
|
|
, m_muSparsity(0.01)
|
|
, m_momentum(0.5)
|
|
, m_doCancel(false)
|
|
, m_useVisibleGaussian(false)
|
|
, m_doRaoBlackwell(false)
|
|
, m_useProbsForHiddenReconstruction(false)
|
|
, m_doSparse(false)
|
|
, m_numGibbs(1)
|
|
{
|
|
Noise_Init(&m_noise, 0x32727155);
|
|
}
|
|
|
|
~Rbm()
|
|
{
|
|
cancel();
|
|
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());
|
|
|
|
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++)
|
|
{
|
|
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();
|
|
|
|
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();
|
|
|
|
} // 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);
|
|
}
|
|
|
|
} // Number of epochs
|
|
}
|
|
|
|
MatrixXd sample(const MatrixXd &src)
|
|
{
|
|
uint32_t i;
|
|
MatrixXd res(src);
|
|
|
|
for (i=0; i < src.array().size(); i++)
|
|
{
|
|
res.array()(i) = (double) src.array()(i) > Noise_Uniform(&m_noise);
|
|
}
|
|
|
|
return res;
|
|
}
|
|
|
|
void train2(const LayerArray<VisibleLayer> &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 vp(m_w.getNumVisible(), batchSize);
|
|
MatrixXd vs(m_w.getNumVisible(), batchSize);
|
|
MatrixXd hp(m_w.getNumHidden(), batchSize);
|
|
MatrixXd hs(m_w.getNumHidden(), batchSize);
|
|
MatrixXd batch(m_w.getNumVisible(), batchSize);
|
|
|
|
VectorXd sumBiasV(m_w.getNumVisible());
|
|
VectorXd sumBiasH(m_w.getNumHidden());
|
|
MatrixXd sumWeights(m_w.getNumVisible(), m_w.getNumHidden());
|
|
|
|
VectorXd deltaBiasV(m_w.getNumVisible());
|
|
VectorXd deltaBiasH(m_w.getNumHidden());
|
|
MatrixXd deltaWeights(m_w.getNumVisible(), m_w.getNumHidden());
|
|
|
|
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.col(i) = vt[t].states();
|
|
}
|
|
|
|
for (epoch=0; epoch < numEpochs; epoch++)
|
|
{
|
|
if (m_doCancel)
|
|
{
|
|
m_doCancel = false;
|
|
break;
|
|
}
|
|
|
|
sumWeights.fill(0);
|
|
sumBiasV.fill(0);
|
|
sumBiasH.fill(0);
|
|
for (i=0; i < batchSize; i++)
|
|
{
|
|
vs = batch.col(i);
|
|
|
|
// h.probsUpdateLogistic(vt[t], m_w, m_lambda, sigma);
|
|
hp = -vs.transpose() * m_w.weights();
|
|
hp.array() = hp.array().exp();
|
|
hp.array() += 1;
|
|
hp.array() = 1.0/hp.array();
|
|
|
|
// Create hidden layer base on training data
|
|
hs = sample(hp);
|
|
|
|
// Update weights (positive phase)
|
|
sumBiasV += vp.colwise().sum();
|
|
if (m_doRaoBlackwell)
|
|
{
|
|
sumWeights += vs * hp.transpose();
|
|
sumBiasH += hp.colwise().sum();
|
|
}
|
|
else
|
|
{
|
|
sumWeights += vs * hs.transpose();
|
|
sumBiasH += hs.colwise().sum();
|
|
}
|
|
|
|
for (gibbs=0; gibbs < m_numGibbs; gibbs++)
|
|
{
|
|
// Create visible reconstruction (a fantasy...)
|
|
if (m_useProbsForHiddenReconstruction)
|
|
{
|
|
vp = -hs * m_w.weights().transpose();
|
|
vp.array() = vp.array().exp();
|
|
vp.array() += 1;
|
|
vp.array() = 1.0/vp.array();
|
|
}
|
|
else
|
|
{
|
|
vs = sample(vp);
|
|
}
|
|
|
|
// Create hidden reconstruction
|
|
hp = -vs.transpose() * m_w.weights();
|
|
hp.array() = hp.array().exp();
|
|
hp.array() += 1;
|
|
hp.array() = 1.0/hp.array();
|
|
}
|
|
|
|
// Update weights (negative phase)
|
|
sumBiasV -= vp.colwise().sum();
|
|
if (m_doRaoBlackwell)
|
|
{
|
|
sumWeights -= vs * hp.transpose();
|
|
sumBiasH -= hp.colwise().sum();
|
|
}
|
|
else
|
|
{
|
|
sumWeights -= vs * hs.transpose();
|
|
sumBiasH -= hs.colwise().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;
|
|
|
|
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);
|
|
}
|
|
|
|
} // Number of epochs
|
|
}
|
|
|
|
double getProgress() const
|
|
{
|
|
return m_progress;
|
|
}
|
|
|
|
double getEnergy(const VectorXd& visible, const VectorXd& hidden)
|
|
{
|
|
double energy;
|
|
|
|
energy = m_w.visibleBias().transpose() * visible;
|
|
energy += m_w.hiddenBias().transpose() * hidden;
|
|
energy += visible.transpose() * m_w.weights() * hidden;
|
|
|
|
return -energy/(m_sigma*m_sigma);
|
|
}
|
|
|
|
void prob(LayerArray<VisibleLayer> &vts)
|
|
{
|
|
uint32_t i, j;
|
|
double z;
|
|
double p;
|
|
|
|
HiddenLayer *h = new HiddenLayer[vts.getSize()];
|
|
|
|
// 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;
|
|
}
|
|
|
|
VectorXd toHidden(const VectorXd& visible)
|
|
{
|
|
HiddenLayer th(m_w.getNumHidden());
|
|
VisibleLayer tv(m_w.getNumVisible(), (const VectorXd*)&visible);
|
|
|
|
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);
|
|
}
|
|
else
|
|
{
|
|
tv.probsUpdateLogistic(th, m_w, m_lambda, m_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();
|
|
}
|
|
|
|
void setSigma(double value)
|
|
{
|
|
m_sigma = value;
|
|
}
|
|
|
|
void setSigmaDecay(double value)
|
|
{
|
|
m_sigmaDecay = value;
|
|
}
|
|
|
|
void setWeightDecay(double value)
|
|
{
|
|
m_weightDecay = value;
|
|
}
|
|
|
|
void setLambda(double value)
|
|
{
|
|
m_lambda = value;
|
|
}
|
|
|
|
void setSparsity(double value)
|
|
{
|
|
m_sparsity = value;
|
|
}
|
|
|
|
void setUseVisibleGaussian(bool flag)
|
|
{
|
|
m_useVisibleGaussian = flag;
|
|
}
|
|
|
|
void setDoRaoBlackwell(bool flag)
|
|
{
|
|
m_doRaoBlackwell = flag;
|
|
}
|
|
|
|
void setUseProbsForHiddenReconstruction(bool flag)
|
|
{
|
|
m_useProbsForHiddenReconstruction = flag;
|
|
}
|
|
|
|
void setDoSparse(bool flag)
|
|
{
|
|
m_doSparse = flag;
|
|
}
|
|
|
|
void setNumGibbs(uint32_t value)
|
|
{
|
|
m_numGibbs = value;
|
|
}
|
|
|
|
void setMuWeights(double value)
|
|
{
|
|
m_muWeights = value;
|
|
}
|
|
|
|
void setMuSparsity(double value)
|
|
{
|
|
m_muSparsity = value;
|
|
}
|
|
|
|
void setMomentum(double value)
|
|
{
|
|
m_momentum = value;
|
|
}
|
|
|
|
void cancel()
|
|
{
|
|
m_doCancel = true;
|
|
// while(m_doCancel);
|
|
}
|
|
|
|
private:
|
|
Weights &m_w;
|
|
RbmListener *m_pListener;
|
|
noise_gen_t m_noise;
|
|
double m_progress;
|
|
double m_sigma;
|
|
double m_sigmaDecay;
|
|
double m_weightDecay;
|
|
double m_lambda;
|
|
double m_sparsity;
|
|
double m_muWeights;
|
|
double m_muSparsity;
|
|
double m_momentum;
|
|
bool m_useVisibleGaussian;
|
|
bool m_doRaoBlackwell;
|
|
bool m_useProbsForHiddenReconstruction;
|
|
bool m_doSparse;
|
|
volatile bool m_doCancel;
|
|
uint32_t m_numGibbs;
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
#endif /* RBM_HPP_ */
|