Want more features on Pastebin? Sign Up, it's FREE!
Guest

Levenberq-Marquardt problem

By: a guest on Mar 4th, 2013  |  syntax: C++  |  size: 8.37 KB  |  views: 89  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include "LevenbergMarquardtTraining.h"
  2.  
  3. #include "NeuralNetwork.h"
  4. #include "ActivationFunctions.h"
  5. #include "../Libraries/alglib/solvers.h"
  6.  
  7. LevenbergMarquardtTraining::LevenbergMarquardtTraining(NeuralNetwork * net)
  8. {
  9.         //efficient algorithm for training neural networks with one hidden layer
  10.         //http://cs.olemiss.edu/~ychen/publications/conference/chen_ijcnn99.pdf
  11.  
  12.         this->net = net;
  13.  
  14.         this->oCount = this->net->GetNeuronsCount(NeuralNetwork::OUTPUT);
  15.         this->hCount = this->net->GetNeuronsCount(NeuralNetwork::HIDDEN);
  16.  
  17.         this->N = this->net->GetWeightsCount();
  18.         this->K = this->oCount;
  19.  
  20.         this->e.resize(this->K);
  21.         this->res.resize(this->N);
  22.        
  23.        
  24.         this->jacobi.setlength(this->K, this->N);
  25.         this->inverse.setlength(this->K, this->K);
  26.  
  27.  
  28.         this->oDelta = Matrix2D(this->oCount, this->oCount);
  29.         this->hDelta = Matrix2D(this->oCount, this->hCount);
  30.  
  31.  
  32. }
  33.  
  34. LevenbergMarquardtTraining::~LevenbergMarquardtTraining()
  35. {
  36. }
  37.  
  38. void LevenbergMarquardtTraining::Train(const std::vector<TrainPair> & train, double errTresshold, int maxIterCount)
  39. {
  40.                
  41.         this->bestWeights.clear();
  42.  
  43.        
  44.  
  45.         double actError = 0;
  46.         double lastError = 0;
  47.         double minError = DBL_MAX;
  48.        
  49.  
  50.         int iterCount = 0;
  51.         double percent = maxIterCount / 100.0;
  52.        
  53.         std::vector<double> weightBckp;
  54.         int m = 1;
  55.         this->mi = 0.01;
  56.  
  57.         //Main loop
  58.         do
  59.         {
  60.                
  61.                 for (uint32 k = 0; k < this->K; k++)
  62.                 {
  63.                         this->e[k] = 0;
  64.                 }
  65.                                                
  66.                
  67.                 //clear jacobi
  68.                 for (uint32 k = 0; k < this->K; k++)
  69.                 {
  70.                         for (uint32 n = 0; n < this->N; n++)
  71.                         {
  72.                                 this->jacobi(k, n) = 0;                
  73.                         }
  74.                 }
  75.                
  76.                 m = 1;
  77.                 actError = 0;
  78.  
  79.                 for (int i = 0; i < train.size(); i++)
  80.                 {                      
  81.                         this->ForwardPass(train[i]);   
  82.                        
  83.                         //calculate deltas for one training input / output
  84.                         this->BackwardPass();
  85.  
  86.  
  87.                         this->CalcJacobi();
  88.  
  89.                         actError += this->TrainError(train[i].output, this->net->GetOutputs());                                                        
  90.                 }
  91.  
  92.                 this->net->GetWeights(weightBckp);
  93.  
  94.                 lastError = actError;
  95.                
  96.  
  97.         JACOBI_UPDATE:
  98.                 actError = 0;
  99.                 this->CalcInverse();
  100.                 this->CalcUpdateRate();
  101.                
  102.                 this->UpdateWeights(weightBckp);
  103.  
  104.                 //calculate error after update
  105.                 //to check, if update was helpfull
  106.                 for (int i = 0; i < train.size(); i++)
  107.                 {                                              
  108.                         this->net->Run(train[i].input);
  109.                         actError += this->TrainError(train[i].output, this->net->GetOutputs());                                                        
  110.                 }
  111.  
  112.                 //according to usefulness of update, do something
  113.                 if (actError <= lastError)
  114.                 {
  115.                         this->mi /= 10;
  116.                 }
  117.                 else if (actError > lastError)
  118.                 {                      
  119.                         if (m <= 5)
  120.                         {
  121.                                 m++;
  122.                                 this->mi *= 10;
  123.                                 this->net->SetWeights(weightBckp);
  124.                                 goto JACOBI_UPDATE;
  125.                         }                      
  126.                 }
  127.                 else if (actError <= errTresshold)
  128.                 {
  129.                         break;
  130.                 }
  131.  
  132.                 if (actError < minError)
  133.                 {                      
  134.                         minError = actError;
  135.                         this->net->GetWeights(this->bestWeights);
  136.                 }
  137.  
  138.                 iterCount++;
  139.  
  140.                 int actPercent = iterCount / percent;
  141.                
  142.                 //if (actPercent % 20 == 0)
  143.                 {
  144.                         //printf("Trained: %d, min-error: %f\n", actPercent, minError);
  145.                         printf("Act-error: %f, min-error: %f\n", actError, minError);
  146.                 }
  147.  
  148.  
  149.         } while (iterCount < maxIterCount);
  150.  
  151.         printf("Error best: %f \n", minError);
  152.        
  153.        
  154.         this->net->GetWeights(this->bestWeights);
  155.  
  156. }
  157.  
  158. void LevenbergMarquardtTraining::ForwardPass(const TrainPair & train)
  159. {
  160.         this->net->Run(train.input);
  161.  
  162.        
  163.         std::vector<double> output = this->net->GetOutputs();
  164.        
  165.        
  166.  
  167.         for (uint32 i = 0; i < this->K; i++)
  168.         {
  169.                 this->e[i] += (train.output[i] - output[i]) * (train.output[i] - output[i]);
  170.         }
  171.        
  172.  
  173. }
  174.  
  175. void LevenbergMarquardtTraining::BackwardPass()
  176. {
  177.         /*
  178.         this->hSlope.clear();
  179.         this->oSlope.clear();
  180.  
  181.         this->net->GetSlope(NeuralNetwork::HIDDEN, this->hSlope);
  182.         this->net->GetSlope(NeuralNetwork::OUTPUT, this->oSlope);
  183.         */
  184.        
  185.         this->oGrads.clear();
  186.         this->hGrads.clear();
  187.         this->oGrads.resize(this->K);
  188.         this->hGrads.resize(this->hCount);
  189.  
  190.         // 1. compute output gradients
  191.         for (uint32 i = 0; i < this->K; i++)
  192.     {                  
  193.                 //this->oGrads[i] = this->oSlope[i];
  194.                 this->oGrads[i] = this->net->hoFunc->DerivateFunction(this->net->hoSums[i]);
  195.     }
  196.        
  197.  
  198.         // 2. compute hidden gradients
  199.         for (uint32 i = 0; i < this->hCount; i++)
  200.     {                          
  201.                 double sum = 0.0;
  202.                
  203.                 for (uint32 j = 0; j < this->K; j++)
  204.                 {
  205.                         sum += this->oGrads[j] * this->net->hoWeights.m[i][j]; // each downstream gradient * outgoing weight
  206.                 }
  207.                 //this->hGrads[i] = sum * this->hSlope[i];             
  208.                 this->hGrads[i] = sum * this->net->ihFunc->DerivateFunction(this->net->ihSums[i]);             
  209.         }
  210.  
  211.  
  212. }
  213.  
  214.  
  215. void LevenbergMarquardtTraining::CalcJacobi()
  216. {
  217.         std::vector<double> output = this->net->GetOutputs();
  218.        
  219.        
  220.        
  221.         //std::vector<double> weights;
  222.         //this->net->GetWeights(weights);
  223.  
  224.         for (uint32 k = 0; k < this->K; k++)
  225.         {
  226.                 uint32 n = 0;
  227.  
  228.                 for (uint32 i = 0; i < this->net->inputCount; i++)
  229.                 {
  230.                         for (uint32 j = 0; j < this->net->hiddenCount; j++)
  231.                         {
  232.                                 //this->jacobi(k, n++) += (/*this->net->inputs[i] * */ output[k] * this->hGrads[j]);
  233.                                 //this->jacobi(k, n) += (this->net->inputs[i] * this->hGrads[j]);
  234.                                 this->jacobi(k, n++) += (this->net->inputs[i] * this->hGrads[j]);
  235.                         }
  236.                 }
  237.  
  238.                 for (uint32 i = 0; i < this->net->hiddenCount; i++)
  239.                 {
  240.                         //this->jacobi(k, n++) += (output[k] * this->hGrads[i]);
  241.                         //this->jacobi(k, n) += (this->hGrads[i]);
  242.                         this->jacobi(k, n++) += (this->hGrads[i]);
  243.                 }
  244.  
  245.  
  246.                 // 4. update hidden to output weights
  247.                 for (uint32 i = 0; i < this->net->hiddenCount; i++)
  248.                 {
  249.                         for (uint32 j = 0; j < this->net->outputCount; j++)
  250.                         {
  251.                                 //this->jacobi(k, n++) += (/*this->net->ihOutputs[i] * */ output[k] * this->oGrads[j]);
  252.                                 //this->jacobi(k, n) += (this->net->ihOutputs[i] * this->oGrads[j]);
  253.                                 this->jacobi(k, n++) += (this->net->ihOutputs[i] * this->oGrads[j]);
  254.                         }
  255.                 }
  256.  
  257.                 // 4b. update hidden to output biases
  258.                 for (uint32 i = 0; i < this->net->outputCount; i++)
  259.                 {
  260.                         //this->jacobi(k, n++) += (output[k] * this->oGrads[i]);       
  261.                         //this->jacobi(k, n) += (this->oGrads[i]);
  262.                         this->jacobi(k, n++) += (this->oGrads[i]);
  263.                 }
  264.  
  265.                 /*
  266.                 for (uint32 n = 0; n < this->N; n++)
  267.                 {
  268.  
  269.                         double grad = 0;
  270.                        
  271.  
  272.  
  273.                         this->jacobi(k, n) += (output[k] * grad);
  274.                         //this->jacobi(k, n) = f->DerivateFunctionFromFx(this->e[k]);
  275.                 }
  276.                 */
  277.         }
  278. }
  279.  
  280.  
  281. void LevenbergMarquardtTraining::CalcInverse()
  282. {
  283.         //fill matrix with values
  284.         //input = I + 1/mi * J * JT
  285.  
  286.         //http://www.alglib.net/translator/man/manual.cpp.html#example_ablas_d_gemm
  287.  
  288.        
  289.         alglib::rmatrixgemm(this->K, this->K, this->N,
  290.                 1.0 / this->mi,
  291.                 this->jacobi, 0, 0, 0,
  292.                 this->jacobi, 0, 0, 1,
  293.                 0,
  294.                 this->inverse, 0, 0);
  295.        
  296.        
  297.         for (uint32 k = 0; k < this->K; k++)
  298.         {
  299.                 this->inverse(k, k) += 1.0; //add I matrix
  300.         }
  301.                
  302.        
  303.         //now compute Inverse
  304.         //http://www.alglib.net/translator/man/manual.cpp.html#sub_rmatrixinverse
  305.         alglib::ae_int_t info;
  306.         alglib::matinvreport inv;      
  307.         alglib::rmatrixinverse(this->inverse, info, inv);
  308.  
  309.        
  310. }
  311.  
  312. void LevenbergMarquardtTraining::CalcUpdateRate()
  313. {
  314.         alglib::real_2d_array tmpNK;
  315.         alglib::real_2d_array tmpNN;
  316.         tmpNK.setlength(this->N, this->K);
  317.         tmpNN.setlength(this->N, this->N);
  318.  
  319.        
  320.  
  321.         // (1/mi * I - 1/mi^2 * JT * Inverse * J) * JT * E
  322.  
  323.         //tmpNK = 1/mi^2 * JT * Inverse
  324.         alglib::rmatrixgemm(this->N, this->K, this->K,
  325.                 1.0 / (this->mi * this->mi),
  326.                 this->jacobi, 0, 0, 1,
  327.                 this->inverse, 0, 0, 0,
  328.                 0,
  329.                 tmpNK, 0, 0);
  330.  
  331.  
  332.         //tmpNN = tmpNK * J
  333.         alglib::rmatrixgemm(this->N, this->N, this->K,
  334.                 1,
  335.                 tmpNK, 0, 0, 0,
  336.                 this->jacobi, 0, 0, 0,
  337.                 0,
  338.                 tmpNN, 0, 0);
  339.  
  340.  
  341.         //tmpNN = 1 / mi * I -  tmpNN
  342.         for (uint32 i = 0; i < this->N; i++)
  343.         {
  344.                 for (uint32 j = 0; j < this->N; j++)
  345.                 {
  346.                         tmpNN(i, j) *= -1;
  347.                 }
  348.         }
  349.        
  350.         for (uint32 n = 0; n < this->N; n++)
  351.         {
  352.                 tmpNN(n, n) += (1.0 / this->mi);
  353.         }
  354.  
  355.  
  356.         //tmpNK = tmpNN * JT
  357.         alglib::rmatrixgemm(this->N, this->K, this->N,
  358.                 1,
  359.                 tmpNN, 0, 0, 0,
  360.                 this->jacobi, 0, 0, 1,
  361.                 0,
  362.                 tmpNK, 0, 0);
  363.  
  364.  
  365.  
  366.         for (uint32 n = 0; n < this->N; n++)
  367.         {
  368.                 this->res[n] = 0;
  369.         }
  370.  
  371.  
  372.         //res = tmpNK * E      
  373.         for (uint32 n = 0; n < this->N; n++)
  374.         {
  375.                 for (uint32 k = 0; k < this->K; k++)
  376.                 {
  377.                         this->res[n] += tmpNK(n, k) * this->e[k];
  378.                 }
  379.         }
  380.  
  381. }
  382.  
  383. void LevenbergMarquardtTraining::UpdateWeights(const std::vector<double> & oldWeights)
  384. {
  385.         std::vector<double> newWeights;
  386.         newWeights.resize(this->N);
  387.  
  388.        
  389.         for (uint32 n = 0; n < this->N; n++)
  390.         {
  391.                 newWeights[n] = oldWeights[n] - this->res[n];          
  392.         }
  393.  
  394.        
  395.         this->mi = 0;
  396.         for (uint32 k = 0; k < this->K; k++)
  397.         {
  398.                 this->mi += this->e[k] * this->e[k];
  399.         }
  400.         this->mi *= 0.01;
  401.        
  402.         this->net->SetWeights(newWeights);
  403. }
clone this paste RAW Paste Data