[mlpack-svn] r10090 - mlpack/trunk/src/mlpack/methods/mog
fastlab-svn at coffeetalk-1.cc.gatech.edu
fastlab-svn at coffeetalk-1.cc.gatech.edu
Mon Oct 31 11:54:32 EDT 2011
Author: jcline3
Date: 2011-10-31 11:54:31 -0400 (Mon, 31 Oct 2011)
New Revision: 10090
Added:
mlpack/trunk/src/mlpack/methods/mog/kmeans.cpp
mlpack/trunk/src/mlpack/methods/mog/kmeans.hpp
mlpack/trunk/src/mlpack/methods/mog/mog_em.cpp
mlpack/trunk/src/mlpack/methods/mog/mog_em.hpp
mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cpp
mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cpp
mlpack/trunk/src/mlpack/methods/mog/mog_l2e.hpp
mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cpp
mlpack/trunk/src/mlpack/methods/mog/optimizers.cpp
mlpack/trunk/src/mlpack/methods/mog/optimizers.hpp
mlpack/trunk/src/mlpack/methods/mog/phi.hpp
Removed:
mlpack/trunk/src/mlpack/methods/mog/kmeans.cc
mlpack/trunk/src/mlpack/methods/mog/kmeans.h
mlpack/trunk/src/mlpack/methods/mog/mog_em.cc
mlpack/trunk/src/mlpack/methods/mog/mog_em.h
mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cc
mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cc
mlpack/trunk/src/mlpack/methods/mog/mog_l2e.h
mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cc
mlpack/trunk/src/mlpack/methods/mog/optimizers.cc
mlpack/trunk/src/mlpack/methods/mog/optimizers.h
mlpack/trunk/src/mlpack/methods/mog/phi.h
Modified:
mlpack/trunk/src/mlpack/methods/mog/CMakeLists.txt
Log:
mog/* to hpp & cpp filenames, fix #defines, etc
Modified: mlpack/trunk/src/mlpack/methods/mog/CMakeLists.txt
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/CMakeLists.txt 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/CMakeLists.txt 2011-10-31 15:54:31 UTC (rev 10090)
@@ -3,18 +3,18 @@
# Define the files we need to compile.
# Anything not in this list will not be compiled into MLPACK.
set(SOURCES
- kmeans.h
- kmeans.cc
+ kmeans.hpp
+ kmeans.cpp
# for mog_em
- mog_em.h
- mog_em.cc
- phi.h
+ mog_em.hpp
+ mog_em.cpp
+ phi.hpp
# for mog_l2e
- mog_l2e.h
- mog_l2e.cc
+ mog_l2e.hpp
+ mog_l2e.cpp
# optimizers
- optimizers.h
- optimizers.cc
+ optimizers.hpp
+ optimizers.cpp
)
# Add directory name to sources.
@@ -28,7 +28,7 @@
# main executable, em
add_executable(mog_em
- mog_em_main.cc
+ mog_em_main.cpp
)
# link dependencies of mog_em
target_link_libraries(mog_em
@@ -37,7 +37,7 @@
# main executable, l2e
add_executable(mog_l2e
- mog_l2e_main.cc
+ mog_l2e_main.cpp
)
target_link_libraries(mog_l2e
mlpack
Deleted: mlpack/trunk/src/mlpack/methods/mog/kmeans.cc
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/kmeans.cc 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/kmeans.cc 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,154 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file kmeans.cc
- *
- * Implementation for the K-means method for getting an initial point.
- *
- */
-#include "kmeans.h"
-
-void KMeans(const arma::mat& data,
- const size_t value_of_k,
- std::vector<arma::vec>& means,
- std::vector<arma::mat>& covars,
- arma::vec& weights) {
- // Set size of vectors and matrices properly.
- means.resize(value_of_k);
- covars.resize(value_of_k);
- for (size_t i = 0; i < value_of_k; i++) {
- means[i].set_size(value_of_k);
- covars[i].set_size(value_of_k, value_of_k);
- }
- weights.set_size(value_of_k);
-
- std::vector<arma::vec> mu, mu_old;
- double* tmpssq = new double[value_of_k];
- double* sig = new double[value_of_k];
- double* sig_best = new double[value_of_k];
- size_t* y = new size_t[value_of_k];
- arma::vec x, diff;
- arma::mat ssq;
- size_t i, j, k, n, t, dim;
- double score, score_old, sum;
-
- n = data.n_cols;
- dim = data.n_rows;
- mu.resize(value_of_k);
- mu_old.resize(value_of_k);
- ssq.set_size(n, value_of_k);
-
- for (i = 0; i < value_of_k; i++) {
- mu[i].set_size(dim);
- mu_old[i].set_size(dim);
- }
-
- x.set_size(dim);
- diff.set_size(dim);
-
- score_old = 999999;
-
- // putting 5 random restarts to obtain the k-means
- for (i = 0; i < 5; i++) {
- t = -1;
- for (k = 0; k < value_of_k; k++){
- t = (t + 1 + (rand() % ((n - 1 - (value_of_k - k)) - (t + 1))));
- mu[k] = data.col(t);
- for(j = 0; j < n; j++) {
- x = data.col(j);
- diff = mu[k] - x;
- ssq(j, k) = dot(diff, diff);
- }
- }
- // This should be an Armadillo function, really.
- double min_val = DBL_MAX;
- for (i = 0; i < ssq.n_rows; i++) {
- for (k = 0; k < ssq.n_cols; k++) {
- if (ssq(i, k) < min_val) {
- min_val = ssq(i, k);
- y[i] = k;
- }
- }
- }
-
- do {
- for (k = 0; k < value_of_k; k++)
- mu_old[k] = mu[k];
-
- for(k = 0; k < value_of_k; k++) {
- size_t p = 0;
- mu[k].zeros();
- for (j = 0; j < n; j++) {
- x = data.col(j);
- if (y[j] == k) {
- mu[k] += x;
- p++;
- }
- }
-
- if (p != 0)
- mu[k] /= p;
-
- for (j = 0; j < n; j++) {
- x = data.col(j);
- diff = mu[k] - x;
- ssq(j, k) = dot(diff, diff);
- }
- }
- // This should be an Armadillo function, really.
- min_val = DBL_MAX;
- for (i = 0; i < ssq.n_rows; i++) {
- for (k = 0; k < ssq.n_cols; k++) {
- if (ssq(i, k) < min_val) {
- min_val = ssq(i, k);
- y[i] = k;
- }
- }
- }
-
- sum = 0;
- for(k = 0; k < value_of_k; k++) {
- diff = mu[k] - mu_old[k];
- sum += dot(diff, diff);
- }
-
- } while (sum != 0);
-
- for (k = 0; k < value_of_k; k++) {
- size_t p = 0;
- tmpssq[k] = 0;
- for (j = 0; j < n; j++) {
- if (y[j] == k) {
- tmpssq[k] += ssq(j, k);
- p++;
- }
- }
- sig[k] = sqrt(tmpssq[k] / p);
- }
-
- score = 0;
- for(k = 0; k < value_of_k; k++) {
- score += tmpssq[k];
- }
- score = score / n;
-
- if (score < score_old) {
- score_old = score;
- for(k = 0; k < value_of_k; k++){
- means[k] = mu[k];
- sig_best[k] = sig[k];
- }
- }
- }
-
- for (k = 0; k < value_of_k; k++) {
- x.fill(sig_best[k]);
- covars[k].diag() = x;
- }
-
- weights.fill(1.0 / value_of_k);
-
- delete[] tmpssq;
- delete[] sig;
- delete[] sig_best;
- delete[] y;
-}
Copied: mlpack/trunk/src/mlpack/methods/mog/kmeans.cpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/kmeans.cc)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/kmeans.cpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/kmeans.cpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,160 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file kmeans.cpp
+ *
+ * Implementation for the K-means method for getting an initial point.
+ *
+ */
+#include "kmeans.hpp"
+
+namespace mlpack {
+namespace gmm {
+
+void KMeans(const arma::mat& data,
+ const size_t value_of_k,
+ std::vector<arma::vec>& means,
+ std::vector<arma::mat>& covars,
+ arma::vec& weights) {
+ // Set size of vectors and matrices properly.
+ means.resize(value_of_k);
+ covars.resize(value_of_k);
+ for (size_t i = 0; i < value_of_k; i++) {
+ means[i].set_size(value_of_k);
+ covars[i].set_size(value_of_k, value_of_k);
+ }
+ weights.set_size(value_of_k);
+
+ std::vector<arma::vec> mu, mu_old;
+ double* tmpssq = new double[value_of_k];
+ double* sig = new double[value_of_k];
+ double* sig_best = new double[value_of_k];
+ size_t* y = new size_t[value_of_k];
+ arma::vec x, diff;
+ arma::mat ssq;
+ size_t i, j, k, n, t, dim;
+ double score, score_old, sum;
+
+ n = data.n_cols;
+ dim = data.n_rows;
+ mu.resize(value_of_k);
+ mu_old.resize(value_of_k);
+ ssq.set_size(n, value_of_k);
+
+ for (i = 0; i < value_of_k; i++) {
+ mu[i].set_size(dim);
+ mu_old[i].set_size(dim);
+ }
+
+ x.set_size(dim);
+ diff.set_size(dim);
+
+ score_old = 999999;
+
+ // putting 5 random restarts to obtain the k-means
+ for (i = 0; i < 5; i++) {
+ t = -1;
+ for (k = 0; k < value_of_k; k++){
+ t = (t + 1 + (rand() % ((n - 1 - (value_of_k - k)) - (t + 1))));
+ mu[k] = data.col(t);
+ for(j = 0; j < n; j++) {
+ x = data.col(j);
+ diff = mu[k] - x;
+ ssq(j, k) = dot(diff, diff);
+ }
+ }
+ // This should be an Armadillo function, really.
+ double min_val = DBL_MAX;
+ for (i = 0; i < ssq.n_rows; i++) {
+ for (k = 0; k < ssq.n_cols; k++) {
+ if (ssq(i, k) < min_val) {
+ min_val = ssq(i, k);
+ y[i] = k;
+ }
+ }
+ }
+
+ do {
+ for (k = 0; k < value_of_k; k++)
+ mu_old[k] = mu[k];
+
+ for(k = 0; k < value_of_k; k++) {
+ size_t p = 0;
+ mu[k].zeros();
+ for (j = 0; j < n; j++) {
+ x = data.col(j);
+ if (y[j] == k) {
+ mu[k] += x;
+ p++;
+ }
+ }
+
+ if (p != 0)
+ mu[k] /= p;
+
+ for (j = 0; j < n; j++) {
+ x = data.col(j);
+ diff = mu[k] - x;
+ ssq(j, k) = dot(diff, diff);
+ }
+ }
+ // This should be an Armadillo function, really.
+ min_val = DBL_MAX;
+ for (i = 0; i < ssq.n_rows; i++) {
+ for (k = 0; k < ssq.n_cols; k++) {
+ if (ssq(i, k) < min_val) {
+ min_val = ssq(i, k);
+ y[i] = k;
+ }
+ }
+ }
+
+ sum = 0;
+ for(k = 0; k < value_of_k; k++) {
+ diff = mu[k] - mu_old[k];
+ sum += dot(diff, diff);
+ }
+
+ } while (sum != 0);
+
+ for (k = 0; k < value_of_k; k++) {
+ size_t p = 0;
+ tmpssq[k] = 0;
+ for (j = 0; j < n; j++) {
+ if (y[j] == k) {
+ tmpssq[k] += ssq(j, k);
+ p++;
+ }
+ }
+ sig[k] = sqrt(tmpssq[k] / p);
+ }
+
+ score = 0;
+ for(k = 0; k < value_of_k; k++) {
+ score += tmpssq[k];
+ }
+ score = score / n;
+
+ if (score < score_old) {
+ score_old = score;
+ for(k = 0; k < value_of_k; k++){
+ means[k] = mu[k];
+ sig_best[k] = sig[k];
+ }
+ }
+ }
+
+ for (k = 0; k < value_of_k; k++) {
+ x.fill(sig_best[k]);
+ covars[k].diag() = x;
+ }
+
+ weights.fill(1.0 / value_of_k);
+
+ delete[] tmpssq;
+ delete[] sig;
+ delete[] sig_best;
+ delete[] y;
+}
+
+}; // namespace gmm
+}; // namespace mlpack
Deleted: mlpack/trunk/src/mlpack/methods/mog/kmeans.h
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/kmeans.h 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/kmeans.h 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,26 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_em.h
- *
- * Defines a Gaussian Mixture model and
- * estimates the parameters of the model
- */
-#ifndef __MLPACK_METHODS_MOG_KMEANS_H
-#define __MLPACK_METHODS_MOG_KMEANS_H
-
-#include <mlpack/core.h>
-
-/**
- * This function computes the k-means of the data and stores the calculated
- * means and covariances in the std::vector of vectors and matrices passed to
- * it. It sets the weights uniformly.
- *
- * This function is used to obtain a starting point for the optimization.
- */
-void KMeans(const arma::mat& data,
- const size_t value_of_k,
- std::vector<arma::vec>& means,
- std::vector<arma::mat>& covars,
- arma::vec& weights);
-
-#endif
Copied: mlpack/trunk/src/mlpack/methods/mog/kmeans.hpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/kmeans.h)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/kmeans.hpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/kmeans.hpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,31 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file kmeans.hpp
+ *
+ * Defines a Gaussian Mixture model and
+ * estimates the parameters of the model
+ */
+#ifndef __MLPACK_METHODS_MOG_KMEANS_HPP
+#define __MLPACK_METHODS_MOG_KMEANS_HPP
+
+#include <mlpack/core.h>
+
+namespace mlpack {
+namespace gmm {
+
+/**
+ * This function computes the k-means of the data and stores the calculated
+ * means and covariances in the std::vector of vectors and matrices passed to
+ * it. It sets the weights uniformly.
+ *
+ * This function is used to obtain a starting point for the optimization.
+ */
+void KMeans(const arma::mat& data,
+ const size_t value_of_k,
+ std::vector<arma::vec>& means,
+ std::vector<arma::mat>& covars,
+ arma::vec& weights);
+}; // namespace gmm
+}; // namespace mlpack
+
+#endif // __MLPACK_METHODS_MOG_KMEANS_HPP
Deleted: mlpack/trunk/src/mlpack/methods/mog/mog_em.cc
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_em.cc 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_em.cc 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,171 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_em.cc
- *
- * Implementation for the loglikelihood function, the EM algorithm
- * and also computes the K-means for getting an initial point
- *
- */
-#include "mog_em.h"
-#include "phi.h"
-#include "kmeans.h"
-
-using namespace mlpack;
-using namespace gmm;
-
-void MoGEM::ExpectationMaximization(const arma::mat& data_points) {
- // Declaration of the variables */
- size_t num_points;
- size_t dim, num_gauss;
- double sum, tmp;
- std::vector<arma::vec> mu_temp, mu;
- std::vector<arma::mat> sigma_temp, sigma;
- arma::vec omega_temp, omega, x;
- arma::mat cond_prob;
- long double l, l_old, best_l, INFTY = 99999, TINY = 1.0e-10;
-
- // Initializing values
- dim = dimension();
- num_gauss = number_of_gaussians();
- num_points = data_points.n_cols;
-
- // Initializing the number of the vectors and matrices
- // according to the parameters input
- mu_temp.resize(num_gauss);
- mu.resize(num_gauss);
- sigma_temp.resize(num_gauss);
- sigma.resize(num_gauss);
- omega_temp.set_size(num_gauss);
- omega.set_size(num_gauss);
-
- // Allocating size to the vectors and matrices
- // according to the dimensionality of the data
- for(size_t i = 0; i < num_gauss; i++) {
- mu_temp[i].set_size(dim);
- mu[i].set_size(dim);
- sigma_temp[i].set_size(dim, dim);
- sigma[i].set_size(dim, dim);
- }
- x.set_size(dim);
- cond_prob.set_size(num_gauss, num_points);
-
- best_l = -INFTY;
- size_t restarts = 0;
- // performing 5 restarts and choosing the best from them
- while (restarts < 5) {
-
- // assign initial values to 'mu', 'sig' and 'omega' using k-means
- KMeans(data_points, num_gauss, mu_temp, sigma_temp, omega_temp);
-
- l_old = -INFTY;
-
- // calculates the loglikelihood value
- l = Loglikelihood(data_points, mu_temp, sigma_temp, omega_temp);
-
- // added a check here to see if any
- // significant change is being made
- // at every iteration
- while (l - l_old > TINY) {
- // calculating the conditional probabilities
- // of choosing a particular gaussian given
- // the data and the present theta value
- for (size_t j = 0; j < num_points; j++) {
- x = data_points.col(j);
- sum = 0;
- for (size_t i = 0; i < num_gauss; i++) {
- tmp = phi(x, mu_temp[i], sigma_temp[i]) * omega_temp[i];
- cond_prob(i, j) = tmp;
- sum += tmp;
- }
- for (size_t i = 0; i < num_gauss; i++) {
- tmp = cond_prob(i, j);
- cond_prob(i, j) = tmp / sum;
- }
- }
-
- // calculating the new value of the mu
- // using the updated conditional probabilities
- for (size_t i = 0; i < num_gauss; i++) {
- sum = 0;
- mu_temp[i].zeros();
- for (size_t j = 0; j < num_points; j++) {
- x = data_points.col(j);
- mu_temp[i] = cond_prob(i, j) * x;
- sum += cond_prob(i, j);
- }
- mu_temp[i] /= sum;
- }
-
- // calculating the new value of the sig
- // using the updated conditional probabilities
- // and the updated mu
- for (size_t i = 0; i < num_gauss; i++) {
- sum = 0;
- sigma_temp[i].zeros();
- for (size_t j = 0; j < num_points; j++) {
- arma::mat co, ro, c;
- c.set_size(dim, dim);
- x = data_points.col(j);
- x -= mu_temp[i];
- c = x * trans(x);
- sigma_temp[i] += cond_prob(i, j) * c;
- sum += cond_prob(i, j);
- }
- sigma_temp[i] /= sum;
- }
-
- // calculating the new values for omega
- // using the updated conditional probabilities
- arma::vec identity_vector;
- identity_vector.set_size(num_points);
- identity_vector = (1.0 / num_points);
- omega_temp = cond_prob * identity_vector;
-
- l_old = l;
- l = Loglikelihood(data_points, mu_temp, sigma_temp, omega_temp);
- }
-
- // putting a check to see if the best one is chosen
- if (l > best_l) {
- best_l = l;
- for (size_t i = 0; i < num_gauss; i++) {
- mu[i] = mu_temp[i];
- sigma[i] = sigma_temp[i];
- }
- omega = omega_temp;
- }
- restarts++;
- }
-
- for (size_t i = 0; i < num_gauss; i++) {
- set_mu(i, mu[i]);
- set_sigma(i, sigma[i]);
- }
- set_omega(omega);
-
- Log::Info << "Log likelihood value of the estimated model: " << best_l << "."
- << std::endl;
- return;
-}
-
-long double MoGEM::Loglikelihood(const arma::mat& data_points,
- const std::vector<arma::vec>& means,
- const std::vector<arma::mat>& covars,
- const arma::vec& weights) {
- size_t i, j;
- arma::vec x;
- long double likelihood, loglikelihood = 0;
-
- x.set_size(data_points.n_rows);
-
- for (j = 0; j < data_points.n_cols; j++) {
- x = data_points.col(j);
- likelihood = 0;
- for(i = 0; i < number_of_gaussians_; i++) {
- likelihood += weights(i) * phi(x, means[i], covars[i]);
- }
- loglikelihood += log(likelihood);
- }
-
- return loglikelihood;
-}
Copied: mlpack/trunk/src/mlpack/methods/mog/mog_em.cpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/mog_em.cc)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_em.cpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_em.cpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,171 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_em.cpp
+ *
+ * Implementation for the loglikelihood function, the EM algorithm
+ * and also computes the K-means for getting an initial point
+ *
+ */
+#include "mog_em.hpp"
+#include "phi.hpp"
+#include "kmeans.hpp"
+
+using namespace mlpack;
+using namespace gmm;
+
+void MoGEM::ExpectationMaximization(const arma::mat& data_points) {
+ // Declaration of the variables */
+ size_t num_points;
+ size_t dim, num_gauss;
+ double sum, tmp;
+ std::vector<arma::vec> mu_temp, mu;
+ std::vector<arma::mat> sigma_temp, sigma;
+ arma::vec omega_temp, omega, x;
+ arma::mat cond_prob;
+ long double l, l_old, best_l, INFTY = 99999, TINY = 1.0e-10;
+
+ // Initializing values
+ dim = dimension();
+ num_gauss = number_of_gaussians();
+ num_points = data_points.n_cols;
+
+ // Initializing the number of the vectors and matrices
+ // according to the parameters input
+ mu_temp.resize(num_gauss);
+ mu.resize(num_gauss);
+ sigma_temp.resize(num_gauss);
+ sigma.resize(num_gauss);
+ omega_temp.set_size(num_gauss);
+ omega.set_size(num_gauss);
+
+ // Allocating size to the vectors and matrices
+ // according to the dimensionality of the data
+ for(size_t i = 0; i < num_gauss; i++) {
+ mu_temp[i].set_size(dim);
+ mu[i].set_size(dim);
+ sigma_temp[i].set_size(dim, dim);
+ sigma[i].set_size(dim, dim);
+ }
+ x.set_size(dim);
+ cond_prob.set_size(num_gauss, num_points);
+
+ best_l = -INFTY;
+ size_t restarts = 0;
+ // performing 5 restarts and choosing the best from them
+ while (restarts < 5) {
+
+ // assign initial values to 'mu', 'sig' and 'omega' using k-means
+ KMeans(data_points, num_gauss, mu_temp, sigma_temp, omega_temp);
+
+ l_old = -INFTY;
+
+ // calculates the loglikelihood value
+ l = Loglikelihood(data_points, mu_temp, sigma_temp, omega_temp);
+
+ // added a check here to see if any
+ // significant change is being made
+ // at every iteration
+ while (l - l_old > TINY) {
+ // calculating the conditional probabilities
+ // of choosing a particular gaussian given
+ // the data and the present theta value
+ for (size_t j = 0; j < num_points; j++) {
+ x = data_points.col(j);
+ sum = 0;
+ for (size_t i = 0; i < num_gauss; i++) {
+ tmp = phi(x, mu_temp[i], sigma_temp[i]) * omega_temp[i];
+ cond_prob(i, j) = tmp;
+ sum += tmp;
+ }
+ for (size_t i = 0; i < num_gauss; i++) {
+ tmp = cond_prob(i, j);
+ cond_prob(i, j) = tmp / sum;
+ }
+ }
+
+ // calculating the new value of the mu
+ // using the updated conditional probabilities
+ for (size_t i = 0; i < num_gauss; i++) {
+ sum = 0;
+ mu_temp[i].zeros();
+ for (size_t j = 0; j < num_points; j++) {
+ x = data_points.col(j);
+ mu_temp[i] = cond_prob(i, j) * x;
+ sum += cond_prob(i, j);
+ }
+ mu_temp[i] /= sum;
+ }
+
+ // calculating the new value of the sig
+ // using the updated conditional probabilities
+ // and the updated mu
+ for (size_t i = 0; i < num_gauss; i++) {
+ sum = 0;
+ sigma_temp[i].zeros();
+ for (size_t j = 0; j < num_points; j++) {
+ arma::mat co, ro, c;
+ c.set_size(dim, dim);
+ x = data_points.col(j);
+ x -= mu_temp[i];
+ c = x * trans(x);
+ sigma_temp[i] += cond_prob(i, j) * c;
+ sum += cond_prob(i, j);
+ }
+ sigma_temp[i] /= sum;
+ }
+
+ // calculating the new values for omega
+ // using the updated conditional probabilities
+ arma::vec identity_vector;
+ identity_vector.set_size(num_points);
+ identity_vector = (1.0 / num_points);
+ omega_temp = cond_prob * identity_vector;
+
+ l_old = l;
+ l = Loglikelihood(data_points, mu_temp, sigma_temp, omega_temp);
+ }
+
+ // putting a check to see if the best one is chosen
+ if (l > best_l) {
+ best_l = l;
+ for (size_t i = 0; i < num_gauss; i++) {
+ mu[i] = mu_temp[i];
+ sigma[i] = sigma_temp[i];
+ }
+ omega = omega_temp;
+ }
+ restarts++;
+ }
+
+ for (size_t i = 0; i < num_gauss; i++) {
+ set_mu(i, mu[i]);
+ set_sigma(i, sigma[i]);
+ }
+ set_omega(omega);
+
+ Log::Info << "Log likelihood value of the estimated model: " << best_l << "."
+ << std::endl;
+ return;
+}
+
+long double MoGEM::Loglikelihood(const arma::mat& data_points,
+ const std::vector<arma::vec>& means,
+ const std::vector<arma::mat>& covars,
+ const arma::vec& weights) {
+ size_t i, j;
+ arma::vec x;
+ long double likelihood, loglikelihood = 0;
+
+ x.set_size(data_points.n_rows);
+
+ for (j = 0; j < data_points.n_cols; j++) {
+ x = data_points.col(j);
+ likelihood = 0;
+ for(i = 0; i < number_of_gaussians_; i++) {
+ likelihood += weights(i) * phi(x, means[i], covars[i]);
+ }
+ loglikelihood += log(likelihood);
+ }
+
+ return loglikelihood;
+}
Deleted: mlpack/trunk/src/mlpack/methods/mog/mog_em.h
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_em.h 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_em.h 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,218 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_em.h
- *
- * Defines a Gaussian Mixture model and
- * estimates the parameters of the model
- */
-#ifndef __MLPACK_METHODS_MOG_MOG_EM_H
-#define __MLPACK_METHODS_MOG_MOG_EM_H
-
-#include <mlpack/core.h>
-
-PARAM_MODULE("mog", "Parameters for the Gaussian mixture model.");
-
-PARAM_INT("k", "The number of Gaussians in the mixture model (defaults to 1).",
- "mog", 1);
-PARAM_INT("d", "The number of dimensions of the data on which the mixture "
- "model is to be fit.", "mog", 0);
-
-namespace mlpack {
-namespace gmm {
-
-/**
- * A Gaussian mixture model class.
- *
- * This class uses maximum likelihood loss functions to
- * estimate the parameters of a gaussian mixture
- * model on a given data via the EM algorithm.
- *
- *
- * Example use:
- *
- * @code
- * MoGEM mog;
- * ArrayList<double> results;
- *
- * mog.Init(number_of_gaussians, dimension);
- * mog.ExpectationMaximization(data, &results, optim_flag);
- * @endcode
- */
-class MoGEM {
- private:
- // The parameters of the mixture model
- std::vector<arma::vec> mu_;
- std::vector<arma::mat> sigma_;
- arma::vec omega_;
- size_t number_of_gaussians_;
- size_t dimension_;
-
- public:
-
- MoGEM() { }
-
- ~MoGEM() { }
-
- void Init(size_t num_gauss, size_t dimension) {
- // Initialize the private variables
- number_of_gaussians_ = num_gauss;
- dimension_ = dimension;
-
- // Resize the ArrayList of Vectors and Matrices
- mu_.resize(number_of_gaussians_);
- sigma_.resize(number_of_gaussians_);
- }
-
- std::vector<arma::vec>& mu() {
- return mu_;
- }
-
- std::vector<arma::mat>& sigma() {
- return sigma_;
- }
-
- arma::vec& omega() {
- return omega_;
- }
-
- size_t number_of_gaussians() {
- return number_of_gaussians_;
- }
-
- size_t dimension() {
- return dimension_;
- }
-
- arma::vec& mu(size_t i) {
- return mu_[i] ;
- }
-
- arma::mat& sigma(size_t i) {
- return sigma_[i];
- }
-
- double omega(size_t i) {
- return omega_[i];
- }
-
- // The set functions
-
- void set_mu(size_t i, arma::vec& mu) {
- assert(i < number_of_gaussians_);
- assert(mu.n_elem == dimension_);
-
- mu_[i] = mu;
- }
-
- void set_sigma(size_t i, arma::mat& sigma) {
- assert(i < number_of_gaussians_);
- assert(sigma.n_rows == dimension_);
- assert(sigma.n_cols == dimension_);
-
- sigma_[i] = sigma;
- }
-
- void set_omega(arma::vec& omega) {
- assert(omega.n_elem == number_of_gaussians());
-
- omega_ = omega;
- return;
- }
-
- /**
- * This function outputs the parameters of the model
- * to an arraylist of doubles
- *
- * @code
- * ArrayList<double> results;
- * mog.OutputResults(&results);
- * @endcode
- */
- void OutputResults(std::vector<double>& results) {
-
- // Initialize the size of the output array
- results.resize(number_of_gaussians_ * (1 + dimension_ * (1 + dimension_)));
-
- // Copy values to the array from the private variables of the class
- for (size_t i = 0; i < number_of_gaussians_; i++) {
- results[i] = omega_[i];
- for (size_t j = 0; j < dimension_; j++) {
- results[number_of_gaussians_ + (i * dimension_) + j] = (mu_[i])[j];
- for (size_t k = 0; k < dimension_; k++) {
- results[number_of_gaussians_ * (1 + dimension_) +
- (i * dimension_ * dimension_) + (j * dimension_) + k] =
- (sigma_[i])(j, k);
- }
- }
- }
- }
-
- /**
- * This function prints the parameters of the model
- *
- * @code
- * mog.Display();
- * @endcode
- */
- void Display() {
- // Output the model parameters as the omega, mu and sigma
- Log::Info << " Omega : [ ";
- for (size_t i = 0; i < number_of_gaussians_; i++) {
- Log::Info << omega_[i] << " ";
- }
- Log::Info << "]" << std::endl;
-
- Log::Info << " Mu : " << std::endl << "[";
- for (size_t i = 0; i < number_of_gaussians_; i++) {
- for (size_t j = 0; j < dimension_ ; j++) {
- Log::Info << (mu_[i])[j];
- }
- Log::Info << ";";
- if (i == (number_of_gaussians_ - 1)) {
- Log::Info << "\b]" << std::endl;
- }
- }
- Log::Info << "Sigma : ";
- for (size_t i = 0; i < number_of_gaussians_; i++) {
- Log::Info << std::endl << "[";
- for (size_t j = 0; j < dimension_ ; j++) {
- for(size_t k = 0; k < dimension_ ; k++) {
- Log::Info << (sigma_[i])(j, k);
- }
- Log::Info << ";";
- }
- Log::Info << "\b]";
- }
- Log::Info << std::endl;
- }
-
- /**
- * This function calculates the parameters of the model
- * using the Maximum Likelihood function via the
- * Expectation Maximization (EM) Algorithm.
- *
- * @code
- * MoG mog;
- * Matrix data = "the data on which you want to fit the model";
- * ArrayList<double> results;
- * mog.ExpectationMaximization(data, &results);
- * @endcode
- */
- void ExpectationMaximization(const arma::mat& data_points);
-
- /**
- * This function computes the loglikelihood of model.
- * This function is used by the 'ExpectationMaximization'
- * function.
- *
- */
- long double Loglikelihood(const arma::mat& data_points,
- const std::vector<arma::vec>& means,
- const std::vector<arma::mat>& covars,
- const arma::vec& weights);
-};
-
-}; // namespace gmm
-}; // namespace mlpack
-
-#endif
Copied: mlpack/trunk/src/mlpack/methods/mog/mog_em.hpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/mog_em.h)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_em.hpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_em.hpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,218 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_em.h
+ *
+ * Defines a Gaussian Mixture model and
+ * estimates the parameters of the model
+ */
+#ifndef __MLPACK_METHODS_MOG_MOG_EM_HPP
+#define __MLPACK_METHODS_MOG_MOG_EM_HPP
+
+#include <mlpack/core.h>
+
+PARAM_MODULE("mog", "Parameters for the Gaussian mixture model.");
+
+PARAM_INT("k", "The number of Gaussians in the mixture model (defaults to 1).",
+ "mog", 1);
+PARAM_INT("d", "The number of dimensions of the data on which the mixture "
+ "model is to be fit.", "mog", 0);
+
+namespace mlpack {
+namespace gmm {
+
+/**
+ * A Gaussian mixture model class.
+ *
+ * This class uses maximum likelihood loss functions to
+ * estimate the parameters of a gaussian mixture
+ * model on a given data via the EM algorithm.
+ *
+ *
+ * Example use:
+ *
+ * @code
+ * MoGEM mog;
+ * ArrayList<double> results;
+ *
+ * mog.Init(number_of_gaussians, dimension);
+ * mog.ExpectationMaximization(data, &results, optim_flag);
+ * @endcode
+ */
+class MoGEM {
+ private:
+ // The parameters of the mixture model
+ std::vector<arma::vec> mu_;
+ std::vector<arma::mat> sigma_;
+ arma::vec omega_;
+ size_t number_of_gaussians_;
+ size_t dimension_;
+
+ public:
+
+ MoGEM() { }
+
+ ~MoGEM() { }
+
+ void Init(size_t num_gauss, size_t dimension) {
+ // Initialize the private variables
+ number_of_gaussians_ = num_gauss;
+ dimension_ = dimension;
+
+ // Resize the ArrayList of Vectors and Matrices
+ mu_.resize(number_of_gaussians_);
+ sigma_.resize(number_of_gaussians_);
+ }
+
+ std::vector<arma::vec>& mu() {
+ return mu_;
+ }
+
+ std::vector<arma::mat>& sigma() {
+ return sigma_;
+ }
+
+ arma::vec& omega() {
+ return omega_;
+ }
+
+ size_t number_of_gaussians() {
+ return number_of_gaussians_;
+ }
+
+ size_t dimension() {
+ return dimension_;
+ }
+
+ arma::vec& mu(size_t i) {
+ return mu_[i] ;
+ }
+
+ arma::mat& sigma(size_t i) {
+ return sigma_[i];
+ }
+
+ double omega(size_t i) {
+ return omega_[i];
+ }
+
+ // The set functions
+
+ void set_mu(size_t i, arma::vec& mu) {
+ assert(i < number_of_gaussians_);
+ assert(mu.n_elem == dimension_);
+
+ mu_[i] = mu;
+ }
+
+ void set_sigma(size_t i, arma::mat& sigma) {
+ assert(i < number_of_gaussians_);
+ assert(sigma.n_rows == dimension_);
+ assert(sigma.n_cols == dimension_);
+
+ sigma_[i] = sigma;
+ }
+
+ void set_omega(arma::vec& omega) {
+ assert(omega.n_elem == number_of_gaussians());
+
+ omega_ = omega;
+ return;
+ }
+
+ /**
+ * This function outputs the parameters of the model
+ * to an arraylist of doubles
+ *
+ * @code
+ * ArrayList<double> results;
+ * mog.OutputResults(&results);
+ * @endcode
+ */
+ void OutputResults(std::vector<double>& results) {
+
+ // Initialize the size of the output array
+ results.resize(number_of_gaussians_ * (1 + dimension_ * (1 + dimension_)));
+
+ // Copy values to the array from the private variables of the class
+ for (size_t i = 0; i < number_of_gaussians_; i++) {
+ results[i] = omega_[i];
+ for (size_t j = 0; j < dimension_; j++) {
+ results[number_of_gaussians_ + (i * dimension_) + j] = (mu_[i])[j];
+ for (size_t k = 0; k < dimension_; k++) {
+ results[number_of_gaussians_ * (1 + dimension_) +
+ (i * dimension_ * dimension_) + (j * dimension_) + k] =
+ (sigma_[i])(j, k);
+ }
+ }
+ }
+ }
+
+ /**
+ * This function prints the parameters of the model
+ *
+ * @code
+ * mog.Display();
+ * @endcode
+ */
+ void Display() {
+ // Output the model parameters as the omega, mu and sigma
+ Log::Info << " Omega : [ ";
+ for (size_t i = 0; i < number_of_gaussians_; i++) {
+ Log::Info << omega_[i] << " ";
+ }
+ Log::Info << "]" << std::endl;
+
+ Log::Info << " Mu : " << std::endl << "[";
+ for (size_t i = 0; i < number_of_gaussians_; i++) {
+ for (size_t j = 0; j < dimension_ ; j++) {
+ Log::Info << (mu_[i])[j];
+ }
+ Log::Info << ";";
+ if (i == (number_of_gaussians_ - 1)) {
+ Log::Info << "\b]" << std::endl;
+ }
+ }
+ Log::Info << "Sigma : ";
+ for (size_t i = 0; i < number_of_gaussians_; i++) {
+ Log::Info << std::endl << "[";
+ for (size_t j = 0; j < dimension_ ; j++) {
+ for(size_t k = 0; k < dimension_ ; k++) {
+ Log::Info << (sigma_[i])(j, k);
+ }
+ Log::Info << ";";
+ }
+ Log::Info << "\b]";
+ }
+ Log::Info << std::endl;
+ }
+
+ /**
+ * This function calculates the parameters of the model
+ * using the Maximum Likelihood function via the
+ * Expectation Maximization (EM) Algorithm.
+ *
+ * @code
+ * MoG mog;
+ * Matrix data = "the data on which you want to fit the model";
+ * ArrayList<double> results;
+ * mog.ExpectationMaximization(data, &results);
+ * @endcode
+ */
+ void ExpectationMaximization(const arma::mat& data_points);
+
+ /**
+ * This function computes the loglikelihood of model.
+ * This function is used by the 'ExpectationMaximization'
+ * function.
+ *
+ */
+ long double Loglikelihood(const arma::mat& data_points,
+ const std::vector<arma::vec>& means,
+ const std::vector<arma::mat>& covars,
+ const arma::vec& weights);
+};
+
+}; // namespace gmm
+}; // namespace mlpack
+
+#endif // __MLPACK_METHODS_MOG_MOG_EM_HPP
Deleted: mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cc
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cc 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cc 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,69 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_em_main.cc
- *
- * This program test drives the parametric estimation
- * of a Gaussian Mixture model using maximum likelihood.
- *
- * PARAMETERS TO BE INPUT:
- *
- * --data
- * This is the file that contains the data on which
- * the model is to be fit
- *
- * --mog_em/K
- * This is the number of gaussians we want to fit
- * on the data, defaults to '1'
- *
- * --output
- * This file will contain the parameters estimated,
- * defaults to 'ouotput.csv'
- *
- */
-#include "mog_em.h"
-
-PROGRAM_INFO("Mixture of Gaussians",
- "This program takes a parametric estimate of a Gaussian mixture model (GMM)"
- " using maximum likelihood.", "mog");
-
-PARAM_STRING_REQ("data", "A file containing the data on which the model has to "
- "be fit.", "mog");
-PARAM_STRING("output", "The file into which the output is to be written into.",
- "mog", "output.csv");
-
-using namespace mlpack;
-using namespace mlpack::gmm;
-
-int main(int argc, char* argv[]) {
- CLI::ParseCommandLine(argc, argv);
-
- ////// READING PARAMETERS AND LOADING DATA //////
- arma::mat data_points;
- data_points.load(CLI::GetParam<std::string>("mog/data").c_str(),
- arma::auto_detect, false, true);
-
- ////// MIXTURE OF GAUSSIANS USING EM //////
- MoGEM mog;
-
- CLI::GetParam<int>("mog/k") = 1;
- CLI::GetParam<int>("mog/d") = data_points.n_rows;
-
- ////// Timing the initialization of the mixture model //////
- CLI::StartTimer("mog/model_init");
- mog.Init(1, data_points.n_rows);
- CLI::StopTimer("mog/model_init");
-
- ////// Computing the parameters of the model using the EM algorithm //////
- std::vector<double> results;
-
- CLI::StartTimer("mog/EM");
- mog.ExpectationMaximization(data_points);
- CLI::StopTimer("mog/EM");
-
- mog.Display();
- mog.OutputResults(results);
-
- ////// OUTPUT RESULTS //////
- // We need a better solution for this. So, currently, we do nothing.
- // XML is probably the right tool for the job.
-}
Copied: mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cc)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_em_main.cpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,69 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_em_main.cpp
+ *
+ * This program test drives the parametric estimation
+ * of a Gaussian Mixture model using maximum likelihood.
+ *
+ * PARAMETERS TO BE INPUT:
+ *
+ * --data
+ * This is the file that contains the data on which
+ * the model is to be fit
+ *
+ * --mog_em/K
+ * This is the number of gaussians we want to fit
+ * on the data, defaults to '1'
+ *
+ * --output
+ * This file will contain the parameters estimated,
+ * defaults to 'ouotput.csv'
+ *
+ */
+#include "mog_em.hpp"
+
+PROGRAM_INFO("Mixture of Gaussians",
+ "This program takes a parametric estimate of a Gaussian mixture model (GMM)"
+ " using maximum likelihood.", "mog");
+
+PARAM_STRING_REQ("data", "A file containing the data on which the model has to "
+ "be fit.", "mog");
+PARAM_STRING("output", "The file into which the output is to be written into.",
+ "mog", "output.csv");
+
+using namespace mlpack;
+using namespace mlpack::gmm;
+
+int main(int argc, char* argv[]) {
+ CLI::ParseCommandLine(argc, argv);
+
+ ////// READING PARAMETERS AND LOADING DATA //////
+ arma::mat data_points;
+ data_points.load(CLI::GetParam<std::string>("mog/data").c_str(),
+ arma::auto_detect, false, true);
+
+ ////// MIXTURE OF GAUSSIANS USING EM //////
+ MoGEM mog;
+
+ CLI::GetParam<int>("mog/k") = 1;
+ CLI::GetParam<int>("mog/d") = data_points.n_rows;
+
+ ////// Timing the initialization of the mixture model //////
+ CLI::StartTimer("mog/model_init");
+ mog.Init(1, data_points.n_rows);
+ CLI::StopTimer("mog/model_init");
+
+ ////// Computing the parameters of the model using the EM algorithm //////
+ std::vector<double> results;
+
+ CLI::StartTimer("mog/EM");
+ mog.ExpectationMaximization(data_points);
+ CLI::StopTimer("mog/EM");
+
+ mog.Display();
+ mog.OutputResults(results);
+
+ ////// OUTPUT RESULTS //////
+ // We need a better solution for this. So, currently, we do nothing.
+ // XML is probably the right tool for the job.
+}
Deleted: mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cc
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cc 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cc 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,328 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_l2e.cc
- *
- * Implementation for L2 loss function, and
- * also some initial points generator
- *
- */
-#include "mog_l2e.h"
-#include "phi.h"
-#include "kmeans.h"
-
-using namespace mlpack;
-using namespace gmm;
-
-long double MoGL2E::L2Error(const arma::mat& data) {
- long double reg, fit, l2e;
-
- reg = RegularizationTerm_();
- fit = GoodnessOfFitTerm_(data);
- l2e = reg - (2 * fit) / data.n_cols;
-
- return l2e;
-}
-
-long double MoGL2E::L2Error(const arma::mat& data, arma::vec& gradients) {
-
- long double reg, fit, l2e;
-
- arma::vec g_reg, g_fit;
- reg = RegularizationTerm_(g_reg);
- fit = GoodnessOfFitTerm_(data, g_fit);
-
- gradients = g_reg - (2 * g_fit) / data.n_cols;
-
- l2e = reg - (2 * fit) / data.n_cols;
-
- return l2e;
-}
-
-long double MoGL2E::RegularizationTerm_() {
- arma::mat phi_mu, sum_covar;
- arma::vec x;
- long double reg, tmpVal;
-
- phi_mu.set_size(number_of_gaussians_, number_of_gaussians_);
- sum_covar.set_size(dimension_, dimension_);
- x = omega_;
-
- for (size_t k = 1; k < number_of_gaussians_; k++) {
- for (size_t j = 0; j < k; j++) {
- sum_covar = sigma_[k] + sigma_[j];
-
- tmpVal = phi(mu_[k], mu_[j], sum_covar);
- phi_mu(j, k) = tmpVal;
- phi_mu(k, j) = tmpVal;
- }
- }
-
- for(size_t k = 0; k < number_of_gaussians_; k++) {
- sum_covar = 2 * sigma_[k];
-
- phi_mu(k, k) = phi(mu_[k], mu_[k], sum_covar);
- }
-
- // Calculating the reg value
- reg = dot(x, x * phi_mu);
-
- return reg;
-}
-
-long double MoGL2E::RegularizationTerm_(arma::vec& g_reg) {
- arma::mat phi_mu, sum_covar;
- arma::vec x, y;
- long double reg, tmpVal;
-
- arma::vec df_dw, g_omega;
- std::vector<arma::vec> g_mu, g_sigma;
- std::vector<std::vector<arma::vec> > dp_d_mu, dp_d_sigma;
-
- phi_mu.set_size(number_of_gaussians_, number_of_gaussians_);
- sum_covar.set_size(dimension_, dimension_);
- x = omega_;
-
- g_mu.resize(number_of_gaussians_);
- g_sigma.resize(number_of_gaussians_);
- dp_d_mu.resize(number_of_gaussians_);
- dp_d_sigma.resize(number_of_gaussians_);
- for(size_t k = 0; k < number_of_gaussians_; k++){
- dp_d_mu[k].resize(number_of_gaussians_);
- dp_d_sigma[k].resize(number_of_gaussians_);
- }
-
- for(size_t k = 1; k < number_of_gaussians_; k++) {
- for(size_t j = 0; j < k; j++) {
- sum_covar = sigma_[k] * sigma_[j];
-
- std::vector<arma::mat> tmp_d_cov;
- arma::vec tmp_dp_d_sigma;
-
- tmp_d_cov.resize(dimension_ * (dimension_ + 1));
-
- for(size_t i = 0; i < (dimension_ * (dimension_ + 1) / 2); i++) {
- tmp_d_cov[i] = (d_sigma_[k])[i];
- tmp_d_cov[(dimension_ * (dimension_ + 1) / 2) + i] = (d_sigma_[j])[i];
- }
-
- tmpVal = phi(mu_[k], mu_[j], sum_covar, tmp_d_cov, dp_d_mu[j][k],
- tmp_dp_d_sigma);
-
- phi_mu(j, k) = tmpVal;
- phi_mu(k, j) = tmpVal;
-
- dp_d_mu[k][j] = -dp_d_mu[j][k];
-
- arma::vec tmp_dp_1(tmp_dp_d_sigma.n_elem / 2);
- arma::vec tmp_dp_2(tmp_dp_d_sigma.n_elem / 2);
- for (size_t i = 0; i < tmp_dp_1.n_elem; i++) {
- tmp_dp_1[i] = tmp_dp_d_sigma[i];
- tmp_dp_2[i] = tmp_dp_d_sigma[(dimension_ * (dimension_ + 1) / 2) + i];
- }
-
- dp_d_sigma[j][k] = tmp_dp_1;
- dp_d_sigma[k][j] = tmp_dp_2;
- }
- }
-
- for (size_t k = 0; k < number_of_gaussians_; k++) {
- sum_covar = 2 * sigma_[k];
-
- arma::vec junk;
- tmpVal = phi(mu_[k], mu_[k], sum_covar, d_sigma_[k], junk,
- dp_d_sigma[k][k]);
-
- phi_mu(k, k) = tmpVal;
-
- dp_d_mu[k][k].zeros(dimension_);
- }
-
- // Calculating the reg value
- reg = dot(x, x * phi_mu);
-
- // Calculating the g_omega values - a vector of size K-1
- df_dw = 2.0 * y;
- g_omega = d_omega_ * df_dw;
-
- // Calculating the g_mu values - K vectors of size D
- for (size_t k = 0; k < number_of_gaussians_; k++) {
- g_mu[k].zeros(dimension_);
-
- for (size_t j = 0; j < number_of_gaussians_; j++) {
- g_mu[k] += x[j] * dp_d_mu[j][k];
- g_mu[k] *= 2.0 * x[k];
- }
-
- // Calculating the g_sigma values - K vectors of size D(D+1)/2
- for (size_t k = 0; k < number_of_gaussians_; k++) {
- g_sigma[k].zeros((dimension_ * (dimension_ + 1)) / 2);
- for (size_t j = 0; j < number_of_gaussians_; j++)
- g_sigma[k] += x[k] * dp_d_sigma[j][k];
- g_sigma[k] *= 2.0 * x[k];
- }
-
- // Making the single gradient vector of size K*(D+1)*(D+2)/2 - 1
- arma::vec tmp_g_reg((number_of_gaussians_ * (dimension_ + 1) *
- (dimension_ * 2) / 2) - 1);
- size_t j = 0;
- for (size_t k = 0; k < g_omega.n_elem; k++)
- tmp_g_reg[k] = g_omega[k];
- j = g_omega.n_elem;
-
- for (size_t k = 0; k < number_of_gaussians_; k++) {
- for (size_t i = 0; i < dimension_; i++)
- tmp_g_reg[j + (k * dimension_) + i] = (g_mu[k])[i];
-
- for(size_t i = 0; i < (dimension_ * (dimension_ + 1) / 2); i++) {
- tmp_g_reg[j + (number_of_gaussians_ * dimension_)
- + k * (dimension_ * (dimension_ + 1) / 2)
- + i] = (g_sigma[k])[i];
- }
- }
-
- g_reg = tmp_g_reg;
- }
-
- return reg;
-}
-
-long double MoGL2E::GoodnessOfFitTerm_(const arma::mat& data) {
- long double fit;
- arma::mat phi_x(number_of_gaussians_, data.n_cols);
- arma::vec identity_vector;
-
- identity_vector.ones(data.n_cols);
-
- for (size_t k = 0; k < number_of_gaussians_; k++)
- for (size_t i = 0; i < data.n_cols; i++)
- phi_x(k, i) = phi(data.unsafe_col(i), mu_[k], sigma_[k]);
-
- fit = dot(omega_ * phi_x, identity_vector);
-
- return fit;
-}
-
-long double MoGL2E::GoodnessOfFitTerm_(const arma::mat& data,
- arma::vec& g_fit) {
- long double fit;
- arma::mat phi_x(number_of_gaussians_, data.n_cols);
- arma::vec weights, x, y, identity_vector;
- arma::vec g_omega,tmp_g_omega;
- std::vector<arma::vec> g_mu, g_sigma;
-
- weights = omega_;
- x.set_size(data.n_rows);
- identity_vector.ones(data.n_cols);
-
- g_mu.resize(number_of_gaussians_);
- g_sigma.resize(number_of_gaussians_);
-
- for(size_t k = 0; k < number_of_gaussians_; k++) {
- g_mu[k].zeros(dimension_);
- g_sigma[k].zeros(dimension_ * (dimension_ + 1) / 2);
-
- for (size_t i = 0; i < data.n_cols; i++) {
- arma::vec tmp_g_mu, tmp_g_sigma;
- phi_x(k, i) = phi(data.unsafe_col(i), mu_[k], sigma_[k], d_sigma_[k],
- tmp_g_mu, tmp_g_sigma);
-
- g_mu[k] += tmp_g_mu;
- g_sigma[k] = tmp_g_sigma;
- }
-
- g_mu[k] *= weights[k];
- g_sigma[k] *= weights[k];
- }
-
- fit = dot(weights * phi_x, identity_vector);
-
- // Calculating the g_omega
- tmp_g_omega = phi_x * identity_vector;
- g_omega = d_omega_ * tmp_g_omega;
-
- // Making the single gradient vector of size K*(D+1)*(D+2)/2
- arma::vec tmp_g_fit((number_of_gaussians_ * (dimension_ + 1) *
- (dimension_ * 2) / 2) - 1);
- size_t j = 0;
- for (size_t k = 0; k < g_omega.n_elem; k++)
- tmp_g_fit[k] = g_omega[k];
- j = g_omega.n_elem;
- for (size_t k = 0; k < number_of_gaussians_; k++) {
- for (size_t i = 0; i < dimension_; i++)
- tmp_g_fit[j + (k * dimension_) + i] = (g_mu[k])[i];
-
- for (size_t i = 0; i < (dimension_ * (dimension_ + 1) / 2); i++)
- tmp_g_fit[j + number_of_gaussians_ * dimension_
- + k * (dimension_ * (dimension_ + 1) / 2) + i] = (g_sigma[k])[i];
- }
-
- g_fit = tmp_g_fit;
-
- return fit;
-}
-
-void MoGL2E::MultiplePointsGenerator(arma::mat& points,
- const arma::mat& d,
- size_t number_of_components) {
-
- size_t i, j, x;
-
- for (i = 0; i < points.n_rows; i++)
- for (j = 0; j < points.n_cols - 1; j++)
- points(i, j) = (rand() % 20001) / 1000 - 10;
-
- for (i = 0; i < points.n_rows; i++) {
- for (j = 0; j < points.n_cols; j++) {
- arma::vec tmp_mu = d.col(rand() % d.n_cols);
- for (x = 0; x < d.n_rows; x++)
- points(i, number_of_components - 1 + (j * d.n_rows) + x) = tmp_mu[x];
- }
- }
-
- for (i = 0; i < points.n_rows; i++)
- for (j = 0; j < points.n_cols; j++)
- for (x = 0; x < (d.n_rows * (d.n_rows + 1) / 2); x++)
- points(i, (number_of_components * (d.n_rows + 1) - 1)
- + (j * (d.n_rows * (d.n_rows + 1) / 2)) + x) = (rand() % 501) / 100;
-
- return;
-}
-
-void MoGL2E::InitialPointGenerator(arma::vec& theta,
- const arma::mat& data,
- size_t k_comp) {
- std::vector<arma::vec> means;
- std::vector<arma::mat> covars;
- arma::vec weights;
- double noise;
-
- weights.set_size(k_comp);
- means.resize(k_comp);
- covars.resize(k_comp);
-
- theta.set_size(k_comp);
-
- for (size_t i = 0; i < k_comp; i++) {
- means[i].set_size(data.n_rows);
- covars[i].set_size(data.n_rows, data.n_rows);
- }
-
- KMeans(data, k_comp, means, covars, weights);
-
- for (size_t k = 0; k < k_comp - 1; k++) {
- noise = (double) (rand() % 10000) / (double) 1000;
- theta[k] = noise - 5;
- }
-
- for (size_t k = 0; k < k_comp; k++) {
- for (size_t j = 0; j < data.n_rows; j++)
- theta[k_comp - 1 + k * data.n_rows + j] = (means[k])[j];
-
- arma::mat u = chol(covars[k]);
- for(size_t j = 0; j < data.n_rows; j++)
- for(size_t i = 0; i < j + 1; i++)
- theta[k_comp - 1 + (k_comp * data.n_rows)
- + (k * data.n_rows * (data.n_rows + 1) / 2)
- + (j * (j + 1) / 2 + i)] = u(i, j) + ((rand() % 501) / 100);
- }
-}
Copied: mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cc)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_l2e.cpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,328 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_l2e.cpp
+ *
+ * Implementation for L2 loss function, and
+ * also some initial points generator
+ *
+ */
+#include "mog_l2e.hpp"
+#include "phi.hpp"
+#include "kmeans.hpp"
+
+using namespace mlpack;
+using namespace gmm;
+
+long double MoGL2E::L2Error(const arma::mat& data) {
+ long double reg, fit, l2e;
+
+ reg = RegularizationTerm_();
+ fit = GoodnessOfFitTerm_(data);
+ l2e = reg - (2 * fit) / data.n_cols;
+
+ return l2e;
+}
+
+long double MoGL2E::L2Error(const arma::mat& data, arma::vec& gradients) {
+
+ long double reg, fit, l2e;
+
+ arma::vec g_reg, g_fit;
+ reg = RegularizationTerm_(g_reg);
+ fit = GoodnessOfFitTerm_(data, g_fit);
+
+ gradients = g_reg - (2 * g_fit) / data.n_cols;
+
+ l2e = reg - (2 * fit) / data.n_cols;
+
+ return l2e;
+}
+
+long double MoGL2E::RegularizationTerm_() {
+ arma::mat phi_mu, sum_covar;
+ arma::vec x;
+ long double reg, tmpVal;
+
+ phi_mu.set_size(number_of_gaussians_, number_of_gaussians_);
+ sum_covar.set_size(dimension_, dimension_);
+ x = omega_;
+
+ for (size_t k = 1; k < number_of_gaussians_; k++) {
+ for (size_t j = 0; j < k; j++) {
+ sum_covar = sigma_[k] + sigma_[j];
+
+ tmpVal = phi(mu_[k], mu_[j], sum_covar);
+ phi_mu(j, k) = tmpVal;
+ phi_mu(k, j) = tmpVal;
+ }
+ }
+
+ for(size_t k = 0; k < number_of_gaussians_; k++) {
+ sum_covar = 2 * sigma_[k];
+
+ phi_mu(k, k) = phi(mu_[k], mu_[k], sum_covar);
+ }
+
+ // Calculating the reg value
+ reg = dot(x, x * phi_mu);
+
+ return reg;
+}
+
+long double MoGL2E::RegularizationTerm_(arma::vec& g_reg) {
+ arma::mat phi_mu, sum_covar;
+ arma::vec x, y;
+ long double reg, tmpVal;
+
+ arma::vec df_dw, g_omega;
+ std::vector<arma::vec> g_mu, g_sigma;
+ std::vector<std::vector<arma::vec> > dp_d_mu, dp_d_sigma;
+
+ phi_mu.set_size(number_of_gaussians_, number_of_gaussians_);
+ sum_covar.set_size(dimension_, dimension_);
+ x = omega_;
+
+ g_mu.resize(number_of_gaussians_);
+ g_sigma.resize(number_of_gaussians_);
+ dp_d_mu.resize(number_of_gaussians_);
+ dp_d_sigma.resize(number_of_gaussians_);
+ for(size_t k = 0; k < number_of_gaussians_; k++){
+ dp_d_mu[k].resize(number_of_gaussians_);
+ dp_d_sigma[k].resize(number_of_gaussians_);
+ }
+
+ for(size_t k = 1; k < number_of_gaussians_; k++) {
+ for(size_t j = 0; j < k; j++) {
+ sum_covar = sigma_[k] * sigma_[j];
+
+ std::vector<arma::mat> tmp_d_cov;
+ arma::vec tmp_dp_d_sigma;
+
+ tmp_d_cov.resize(dimension_ * (dimension_ + 1));
+
+ for(size_t i = 0; i < (dimension_ * (dimension_ + 1) / 2); i++) {
+ tmp_d_cov[i] = (d_sigma_[k])[i];
+ tmp_d_cov[(dimension_ * (dimension_ + 1) / 2) + i] = (d_sigma_[j])[i];
+ }
+
+ tmpVal = phi(mu_[k], mu_[j], sum_covar, tmp_d_cov, dp_d_mu[j][k],
+ tmp_dp_d_sigma);
+
+ phi_mu(j, k) = tmpVal;
+ phi_mu(k, j) = tmpVal;
+
+ dp_d_mu[k][j] = -dp_d_mu[j][k];
+
+ arma::vec tmp_dp_1(tmp_dp_d_sigma.n_elem / 2);
+ arma::vec tmp_dp_2(tmp_dp_d_sigma.n_elem / 2);
+ for (size_t i = 0; i < tmp_dp_1.n_elem; i++) {
+ tmp_dp_1[i] = tmp_dp_d_sigma[i];
+ tmp_dp_2[i] = tmp_dp_d_sigma[(dimension_ * (dimension_ + 1) / 2) + i];
+ }
+
+ dp_d_sigma[j][k] = tmp_dp_1;
+ dp_d_sigma[k][j] = tmp_dp_2;
+ }
+ }
+
+ for (size_t k = 0; k < number_of_gaussians_; k++) {
+ sum_covar = 2 * sigma_[k];
+
+ arma::vec junk;
+ tmpVal = phi(mu_[k], mu_[k], sum_covar, d_sigma_[k], junk,
+ dp_d_sigma[k][k]);
+
+ phi_mu(k, k) = tmpVal;
+
+ dp_d_mu[k][k].zeros(dimension_);
+ }
+
+ // Calculating the reg value
+ reg = dot(x, x * phi_mu);
+
+ // Calculating the g_omega values - a vector of size K-1
+ df_dw = 2.0 * y;
+ g_omega = d_omega_ * df_dw;
+
+ // Calculating the g_mu values - K vectors of size D
+ for (size_t k = 0; k < number_of_gaussians_; k++) {
+ g_mu[k].zeros(dimension_);
+
+ for (size_t j = 0; j < number_of_gaussians_; j++) {
+ g_mu[k] += x[j] * dp_d_mu[j][k];
+ g_mu[k] *= 2.0 * x[k];
+ }
+
+ // Calculating the g_sigma values - K vectors of size D(D+1)/2
+ for (size_t k = 0; k < number_of_gaussians_; k++) {
+ g_sigma[k].zeros((dimension_ * (dimension_ + 1)) / 2);
+ for (size_t j = 0; j < number_of_gaussians_; j++)
+ g_sigma[k] += x[k] * dp_d_sigma[j][k];
+ g_sigma[k] *= 2.0 * x[k];
+ }
+
+ // Making the single gradient vector of size K*(D+1)*(D+2)/2 - 1
+ arma::vec tmp_g_reg((number_of_gaussians_ * (dimension_ + 1) *
+ (dimension_ * 2) / 2) - 1);
+ size_t j = 0;
+ for (size_t k = 0; k < g_omega.n_elem; k++)
+ tmp_g_reg[k] = g_omega[k];
+ j = g_omega.n_elem;
+
+ for (size_t k = 0; k < number_of_gaussians_; k++) {
+ for (size_t i = 0; i < dimension_; i++)
+ tmp_g_reg[j + (k * dimension_) + i] = (g_mu[k])[i];
+
+ for(size_t i = 0; i < (dimension_ * (dimension_ + 1) / 2); i++) {
+ tmp_g_reg[j + (number_of_gaussians_ * dimension_)
+ + k * (dimension_ * (dimension_ + 1) / 2)
+ + i] = (g_sigma[k])[i];
+ }
+ }
+
+ g_reg = tmp_g_reg;
+ }
+
+ return reg;
+}
+
+long double MoGL2E::GoodnessOfFitTerm_(const arma::mat& data) {
+ long double fit;
+ arma::mat phi_x(number_of_gaussians_, data.n_cols);
+ arma::vec identity_vector;
+
+ identity_vector.ones(data.n_cols);
+
+ for (size_t k = 0; k < number_of_gaussians_; k++)
+ for (size_t i = 0; i < data.n_cols; i++)
+ phi_x(k, i) = phi(data.unsafe_col(i), mu_[k], sigma_[k]);
+
+ fit = dot(omega_ * phi_x, identity_vector);
+
+ return fit;
+}
+
+long double MoGL2E::GoodnessOfFitTerm_(const arma::mat& data,
+ arma::vec& g_fit) {
+ long double fit;
+ arma::mat phi_x(number_of_gaussians_, data.n_cols);
+ arma::vec weights, x, y, identity_vector;
+ arma::vec g_omega,tmp_g_omega;
+ std::vector<arma::vec> g_mu, g_sigma;
+
+ weights = omega_;
+ x.set_size(data.n_rows);
+ identity_vector.ones(data.n_cols);
+
+ g_mu.resize(number_of_gaussians_);
+ g_sigma.resize(number_of_gaussians_);
+
+ for(size_t k = 0; k < number_of_gaussians_; k++) {
+ g_mu[k].zeros(dimension_);
+ g_sigma[k].zeros(dimension_ * (dimension_ + 1) / 2);
+
+ for (size_t i = 0; i < data.n_cols; i++) {
+ arma::vec tmp_g_mu, tmp_g_sigma;
+ phi_x(k, i) = phi(data.unsafe_col(i), mu_[k], sigma_[k], d_sigma_[k],
+ tmp_g_mu, tmp_g_sigma);
+
+ g_mu[k] += tmp_g_mu;
+ g_sigma[k] = tmp_g_sigma;
+ }
+
+ g_mu[k] *= weights[k];
+ g_sigma[k] *= weights[k];
+ }
+
+ fit = dot(weights * phi_x, identity_vector);
+
+ // Calculating the g_omega
+ tmp_g_omega = phi_x * identity_vector;
+ g_omega = d_omega_ * tmp_g_omega;
+
+ // Making the single gradient vector of size K*(D+1)*(D+2)/2
+ arma::vec tmp_g_fit((number_of_gaussians_ * (dimension_ + 1) *
+ (dimension_ * 2) / 2) - 1);
+ size_t j = 0;
+ for (size_t k = 0; k < g_omega.n_elem; k++)
+ tmp_g_fit[k] = g_omega[k];
+ j = g_omega.n_elem;
+ for (size_t k = 0; k < number_of_gaussians_; k++) {
+ for (size_t i = 0; i < dimension_; i++)
+ tmp_g_fit[j + (k * dimension_) + i] = (g_mu[k])[i];
+
+ for (size_t i = 0; i < (dimension_ * (dimension_ + 1) / 2); i++)
+ tmp_g_fit[j + number_of_gaussians_ * dimension_
+ + k * (dimension_ * (dimension_ + 1) / 2) + i] = (g_sigma[k])[i];
+ }
+
+ g_fit = tmp_g_fit;
+
+ return fit;
+}
+
+void MoGL2E::MultiplePointsGenerator(arma::mat& points,
+ const arma::mat& d,
+ size_t number_of_components) {
+
+ size_t i, j, x;
+
+ for (i = 0; i < points.n_rows; i++)
+ for (j = 0; j < points.n_cols - 1; j++)
+ points(i, j) = (rand() % 20001) / 1000 - 10;
+
+ for (i = 0; i < points.n_rows; i++) {
+ for (j = 0; j < points.n_cols; j++) {
+ arma::vec tmp_mu = d.col(rand() % d.n_cols);
+ for (x = 0; x < d.n_rows; x++)
+ points(i, number_of_components - 1 + (j * d.n_rows) + x) = tmp_mu[x];
+ }
+ }
+
+ for (i = 0; i < points.n_rows; i++)
+ for (j = 0; j < points.n_cols; j++)
+ for (x = 0; x < (d.n_rows * (d.n_rows + 1) / 2); x++)
+ points(i, (number_of_components * (d.n_rows + 1) - 1)
+ + (j * (d.n_rows * (d.n_rows + 1) / 2)) + x) = (rand() % 501) / 100;
+
+ return;
+}
+
+void MoGL2E::InitialPointGenerator(arma::vec& theta,
+ const arma::mat& data,
+ size_t k_comp) {
+ std::vector<arma::vec> means;
+ std::vector<arma::mat> covars;
+ arma::vec weights;
+ double noise;
+
+ weights.set_size(k_comp);
+ means.resize(k_comp);
+ covars.resize(k_comp);
+
+ theta.set_size(k_comp);
+
+ for (size_t i = 0; i < k_comp; i++) {
+ means[i].set_size(data.n_rows);
+ covars[i].set_size(data.n_rows, data.n_rows);
+ }
+
+ KMeans(data, k_comp, means, covars, weights);
+
+ for (size_t k = 0; k < k_comp - 1; k++) {
+ noise = (double) (rand() % 10000) / (double) 1000;
+ theta[k] = noise - 5;
+ }
+
+ for (size_t k = 0; k < k_comp; k++) {
+ for (size_t j = 0; j < data.n_rows; j++)
+ theta[k_comp - 1 + k * data.n_rows + j] = (means[k])[j];
+
+ arma::mat u = chol(covars[k]);
+ for(size_t j = 0; j < data.n_rows; j++)
+ for(size_t i = 0; i < j + 1; i++)
+ theta[k_comp - 1 + (k_comp * data.n_rows)
+ + (k * data.n_rows * (data.n_rows + 1) / 2)
+ + (j * (j + 1) / 2 + i)] = u(i, j) + ((rand() % 501) / 100);
+ }
+}
Deleted: mlpack/trunk/src/mlpack/methods/mog/mog_l2e.h
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_l2e.h 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_l2e.h 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,504 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_l2e.h
- *
- * Defines a Gaussian Mixture model and
- * estimates the parameters of the model
- *
- */
-#ifndef __MLPACK_METHODS_MOG_MOG_L2E_H
-#define __MLPACK_METHODS_MOG_MOG_L2E_H
-
-#include <mlpack/core.h>
-
-PARAM_MODULE("mog_l2e", "Parameters for the Gaussian mixture model.");
-
-PARAM_INT("k", "The number of Gaussians in the mixture model (defaults to 1).",
- "mog_l2e", 1);
-PARAM_INT("d", "The number of dimensions of the data on which the mixture "
- "model is to be fit.", "mog_l2e", 0);
-
-namespace mlpack {
-namespace gmm {
-
-/**
- * A Gaussian mixture model class.
- *
- * This class uses L2 loss function to
- * estimate the parameters of a gaussian mixture
- * model on a given data.
- *
- * The parameters are converted for optimization
- * to maintain the following facts:
- * - the weights sum to one
- * - for this, the weights were parameterized using
- * the logistic function
- * - the covariance matrix is always positive definite
- * - for this, the Cholesky decomposition is used
- *
- * Example use:
- *
- * @code
- * MoGL2E mog;
- * ArrayList<double> results;
- * double *params;
- *
- * mog.MakeModel(number_of_gaussians, dimension, params);
- * mog.L2Error(data);
- * mog.OutputResults(&results);
- * @endcode
- */
-class MoGL2E {
- private:
- // The parameters of the Mixture model
- std::vector<arma::vec> mu_;
- std::vector<arma::mat> sigma_;
- arma::vec omega_;
- size_t number_of_gaussians_;
- size_t dimension_;
-
- // The differential for the paramterization
- // for optimization
- arma::mat d_omega_;
- std::vector<std::vector<arma::mat> > d_sigma_;
-
- public:
-
- MoGL2E() { }
-
- ~MoGL2E() { }
-
- void Init(size_t num_gauss, size_t dimension) {
- // Destruct everything to initialize afresh
- mu_.clear();
- sigma_.clear();
- d_sigma_.clear();
-
- // Initialize the private variables
- number_of_gaussians_ = num_gauss;
- dimension_ = dimension;
-
- // Resize the vector of vectors and matrices
- mu_.resize(number_of_gaussians_);
- sigma_.resize(number_of_gaussians_);
- }
-
- void Resize_d_sigma_() {
- d_sigma_.resize(number_of_gaussians_);
- for(size_t i = 0; i < number_of_gaussians_; i++)
- d_sigma_[i].resize(dimension_ * (dimension_ + 1) / 2);
- }
-
- /**
- * This function uses the parameters used for optimization
- * and converts it into athe parameters of a Gaussian
- * mixture model. This is to be used when you do not want
- * the gradient values.
- *
- * Example use:
- *
- * @code
- * MoGL2E mog;
- * mog.MakeModel(number_of_gaussians, dimension,
- * parameters_for_optimization);
- * @endcode
- */
-
- void MakeModel(size_t num_mods, size_t dimension, const arma::vec& theta) {
- arma::vec temp_mu(dimension);
- arma::mat lower_triangle_matrix;
- double sum, s_min = 0.01;
-
- Init(num_mods, dimension);
- lower_triangle_matrix.set_size(dimension, dimension);
-
- // calculating the omega values
- arma::vec temp_array = exp(theta);
- temp_array[num_mods - 1] = 1;
- sum = accu(temp_array);
-
- temp_array /= sum;
- set_omega(temp_array);
-
- // calculating the mu values
- for(size_t k = 0; k < num_mods; k++) {
- for(size_t j = 0; j < dimension; j++) {
- temp_mu[j] = theta[num_mods + (k * dimension) + j - 1];
- }
- set_mu(k, temp_mu);
- }
-
- // calculating the sigma values
- // using a lower triangular matrix and its transpose
- // to obtain a positive definite symmetric matrix
- arma::mat sigma_temp(dimension, dimension);
- for(size_t k = 0; k < num_mods; k++) {
- lower_triangle_matrix.zeros();
- for(size_t j = 0; j < dimension; j++) {
- for(size_t i = 0; i < j; i++) {
- lower_triangle_matrix(j, i) = theta[(num_mods - 1)
- + (num_mods * dimension) + k * (dimension * (dimension + 1) / 2)
- + (j * (j + 1) / 2) + i];
- }
- lower_triangle_matrix(j, j) = theta[(num_mods - 1)
- + (num_mods * dimension) + k * (dimension * (dimension + 1) / 2)
- + (j * (j + 1) / 2) + j] + s_min;
- }
- sigma_temp = lower_triangle_matrix * trans(lower_triangle_matrix);
- set_sigma(k, sigma_temp);
- }
- }
-
- /**
- * This function uses the parameters used for optimization
- * and converts it into the parameters of a Gaussian
- * mixture model. This is to be used when you want
- * the gradient values.
- *
- * Example use:
- *
- * @code
- * MoGL2E mog;
- * mog.MakeModelWithGradients(number_of_gaussians, dimension,
- * parameters_for_optimization);
- * @endcode
- */
- void MakeModelWithGradients(size_t num_mods,
- size_t dimension,
- const arma::vec& theta) {
- arma::vec temp_mu(dimension);
- arma::mat lower_triangle_matrix;
- double sum, s_min = 0.01;
-
- Init(num_mods, dimension);
- lower_triangle_matrix.set_size(dimension, dimension);
-
- // calculating the omega values
- arma::vec temp_array = exp(theta);
- temp_array[num_mods - 1] = 1;
- sum = accu(temp_array);
-
- temp_array /= sum;
- set_omega(temp_array);
-
- // calculating the d_omega values
- arma::mat d_omega_temp(num_mods - 1, num_mods);
- d_omega_temp.zeros();
- for (size_t i = 0; i < num_mods - 1; i++) {
- for (size_t j = 0; j < i; j++) {
- d_omega_temp(i, j) = -(omega_[i] * omega_[j]);
- d_omega_temp(j, i) = -(omega_[i] * omega_[j]);
- }
- d_omega_temp(i, i) = omega_[i] * (1 - omega_[i]);
- }
-
- for (size_t i = 0; i < num_mods - 1; i++)
- d_omega_temp(i, num_mods - 1) = -(omega_[i] * omega_[num_mods - 1]);
-
- set_d_omega(d_omega_temp);
-
- // calculating the mu values
- for (size_t k = 0; k < num_mods; k++) {
- for (size_t j = 0; j < dimension; j++)
- temp_mu[j] = theta[num_mods + (k * dimension) + j - 1];
- set_mu(k, temp_mu);
- }
- // d_mu is not computed because it is implicitly known
- // since no parameterization is applied on them
-
- // using a lower triangular matrix and its transpose
- // to obtain a positive definite symmetric matrix
-
- // initializing the d_sigma values
-
- arma::mat d_sigma_temp(dimension, dimension);
- Resize_d_sigma_();
-
- // calculating the sigma values
- arma::mat sigma_temp(dimension, dimension);
- for (size_t k = 0; k < num_mods; k++) {
- lower_triangle_matrix.zeros();
- for (size_t j = 0; j < dimension; j++) {
- for (size_t i = 0; i < j; i++) {
- lower_triangle_matrix(j, i) = theta[(num_mods - 1)
- + (num_mods * dimension) + k * (dimension * (dimension + 1) / 2)
- + (j * (j + 1) / 2) + i];
- }
- lower_triangle_matrix(j, j) = theta[(num_mods - 1)
- + (num_mods * dimension) + k * (dimension * (dimension + 1) / 2)
- + (j * (j + 1) / 2) + j] + s_min;
- }
- sigma_temp = lower_triangle_matrix * trans(lower_triangle_matrix);
- set_sigma(k, sigma_temp);
-
- // calculating the d_sigma values
- for (size_t i = 0; i < dimension; i++) {
- for (size_t in = 0; in < i + 1; in++) {
- arma::mat d_sigma_d_r(dimension, dimension);
- d_sigma_d_r.zeros();
- d_sigma_d_r(i, in) = 1.0;
-
- d_sigma_temp = d_sigma_d_r * trans(lower_triangle_matrix) +
- lower_triangle_matrix * trans(d_sigma_d_r);
- set_d_sigma(k, (i * (i + 1) / 2) + in, d_sigma_temp);
- }
- }
- }
- }
-
- ////// THE GET FUNCTCLINS //////
- std::vector<arma::vec>& mu() {
- return mu_;
- }
-
- std::vector<arma::mat>& sigma() {
- return sigma_;
- }
-
- arma::vec& omega() {
- return omega_;
- }
-
- size_t number_of_gaussians() {
- return number_of_gaussians_;
- }
-
- size_t dimension() {
- return dimension_;
- }
-
- arma::vec& mu(size_t i) {
- return mu_[i] ;
- }
-
- arma::mat& sigma(size_t i) {
- return sigma_[i];
- }
-
- double omega(size_t i) {
- return omega_[i];
- }
-
- arma::mat& d_omega() {
- return d_omega_;
- }
-
- std::vector<std::vector<arma::mat> >& d_sigma() {
- return d_sigma_;
- }
-
- std::vector<arma::mat>& d_sigma(size_t i) {
- return d_sigma_[i];
- }
-
- ////// THE SET FUNCTCLINS //////
-
- void set_mu(size_t i, const arma::vec& mu) {
- assert(i < number_of_gaussians_);
- assert(mu.n_elem == dimension_);
-
- mu_[i] = mu;
- }
-
- void set_sigma(size_t i, const arma::mat& sigma) {
- assert(i < number_of_gaussians_);
- assert(sigma.n_rows == dimension_);
- assert(sigma.n_cols == dimension_);
-
- sigma_[i] = sigma;
- }
-
- void set_omega(const arma::vec& omega) {
- assert(omega.n_elem == number_of_gaussians_);
- omega_ = omega;
- }
-
- void set_d_omega(const arma::mat& d_omega) {
- d_omega_ = d_omega;
- }
-
- void set_d_sigma(size_t i, size_t j, const arma::mat& d_sigma_i_j) {
- d_sigma_[i][j] = d_sigma_i_j;
- }
-
- /**
- * This function outputs the parameters of the model
- * to an arraylist of doubles
- *
- * @code
- * ArrayList<double> results;
- * mog.OutputResults(&results);
- * @endcode
- */
- void OutputResults(std::vector<double>& results) {
-
- // Initialize the size of the output array
- results.resize(number_of_gaussians_ * (1 + dimension_ * (1 + dimension_)));
-
- // Copy values to the array from the private variables of the class
- for (size_t i = 0; i < number_of_gaussians_; i++) {
- results[i] = omega_[i];
- for (size_t j = 0; j < dimension_; j++) {
- results[number_of_gaussians_ + (i * dimension_) + j] = (mu_[i])[j];
- for (size_t k = 0; k < dimension_; k++) {
- results[number_of_gaussians_ * (1 + dimension_)
- + (i * dimension_ * dimension_) + (j * dimension_)
- + k] = (sigma_[i])(j, k);
- }
- }
- }
- }
-
- /**
- * This function prints the parameters of the model
- *
- * @code
- * mog.Display();
- * @endcode
- */
- void Display() {
- // Output the model parameters as the omega, mu and sigma
- Log::Info << " Omega : [ ";
- for (size_t i = 0; i < number_of_gaussians_; i++) {
- Log::Info << omega_[i];
- }
- Log::Info << "]" << std::endl;
- Log::Info << " Mu : " << std::endl << "[";
- for (size_t i = 0; i < number_of_gaussians_; i++) {
- for (size_t j = 0; j < dimension_ ; j++) {
- Log::Info << (mu_[i])[j];
- }
- Log::Info << ";";
- if (i == (number_of_gaussians_ - 1)) {
- Log::Info << "\b]" << std::endl;
- }
- }
- Log::Info << "Sigma : ";
- for (size_t i = 0; i < number_of_gaussians_; i++) {
- Log::Info << std::endl << "[";
- for (size_t j = 0; j < dimension_ ; j++) {
- for (size_t k = 0; k < dimension_ ; k++) {
- Log::Info << (sigma_[i])(j, k);
- }
- Log::Info << ";";
- }
- Log::Info << "\b]";
- }
- Log::Info << std::endl;
- }
-
- /**
- * This function calculates the L2 error and
- * the gradient of the error with respect to the
- * parameters given the data and the parameterized
- * mixture
- *
- * Example use:
- *
- * @code
- * const Matrix data;
- * MoGL2E mog;
- * size_t num_gauss, dimension;
- * double *params; // get the parameters
- *
- * mog.MakeModel(num_gauss, dimension, params);
- * mog.L2Error(data);
- * @endcode
- */
- long double L2Error(const arma::mat& data);
- long double L2Error(const arma::mat& data, arma::vec& gradients);
-
- /**
- * Calculates the regularization value for a
- * Gaussian mixture and its gradient with
- * respect to the parameters
- *
- * Used by the 'L2Error' function to calculate
- * the regularization part of the error
- */
- long double RegularizationTerm_();
- long double RegularizationTerm_(arma::vec& g_reg);
-
- /**
- * Calculates the goodness-of-fit value for a
- * Gaussian mixture and its gradient with
- * respect to the parameters
- *
- * Used by the 'L2Error' function to calculate
- * the goodness-of-fit part of the error
- */
- long double GoodnessOfFitTerm_(const arma::mat& data);
- long double GoodnessOfFitTerm_(const arma::mat& data, arma::vec& g_fit);
-
- /**
- * This function computes multiple number of starting points
- * required for the Nelder Mead method
- *
- * Example use:
- * @code
- * double **p;
- * size_t n, num_gauss;
- * const Matrix data;
- *
- * MoGL2E::MultiplePointsGeneratot(p, n, data, num_gauss);
- * @endcode
- */
- static void MultiplePointsGenerator(arma::mat& points,
- const arma::mat& d,
- size_t number_of_components);
-
- /**
- * This function parameterizes the starting point obtained
- * from the 'k_means" for optimization purposes using the
- * Quasi Newton method
- *
- * Example use:
- * @code
- * double *p;
- * size_t num_gauss;
- * const Matrix data;
- *
- * MoGL2E::InitialPointGeneratot(p, data, num_gauss);
- * @endcode
- */
- static void InitialPointGenerator(arma::vec& theta,
- const arma::mat& data,
- size_t k_comp);
-
- /**
- * This is the function which would be used for
- * optimization. It creates its own object of
- * class MoGL2E and returns the L2 error
- * and the gradient which are computed by
- * the functions of the class
- *
- */
- static long double L2ErrorForOpt(const arma::vec& params,
- const arma::mat& data) {
- MoGL2E model;
- size_t num_gauss = (params.n_elem + 1) * 2 /
- ((data.n_rows + 1) * (data.n_rows + 2));
-
- model.MakeModel(num_gauss, data.n_rows, params);
-
- return model.L2Error(data);
- }
-
- static long double L2ErrorForOpt(const arma::vec& params,
- const arma::mat& data,
- arma::vec& gradient) {
-
- MoGL2E model;
- size_t num_gauss = (params.n_elem + 1) * 2 /
- ((data.n_rows + 1) * (data.n_rows + 2));
-
- model.MakeModelWithGradients(num_gauss, data.n_rows, params);
-
- return model.L2Error(data, gradient);
- }
-};
-
-}; // namespace gmm
-}; // namespace mlpack
-
-#endif
Copied: mlpack/trunk/src/mlpack/methods/mog/mog_l2e.hpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/mog_l2e.h)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_l2e.hpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_l2e.hpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,504 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_l2e.hpp
+ *
+ * Defines a Gaussian Mixture model and
+ * estimates the parameters of the model
+ *
+ */
+#ifndef __MLPACK_METHODS_MOG_MOG_L2E_HPP
+#define __MLPACK_METHODS_MOG_MOG_L2E_HPP
+
+#include <mlpack/core.h>
+
+PARAM_MODULE("mog_l2e", "Parameters for the Gaussian mixture model.");
+
+PARAM_INT("k", "The number of Gaussians in the mixture model (defaults to 1).",
+ "mog_l2e", 1);
+PARAM_INT("d", "The number of dimensions of the data on which the mixture "
+ "model is to be fit.", "mog_l2e", 0);
+
+namespace mlpack {
+namespace gmm {
+
+/**
+ * A Gaussian mixture model class.
+ *
+ * This class uses L2 loss function to
+ * estimate the parameters of a gaussian mixture
+ * model on a given data.
+ *
+ * The parameters are converted for optimization
+ * to maintain the following facts:
+ * - the weights sum to one
+ * - for this, the weights were parameterized using
+ * the logistic function
+ * - the covariance matrix is always positive definite
+ * - for this, the Cholesky decomposition is used
+ *
+ * Example use:
+ *
+ * @code
+ * MoGL2E mog;
+ * ArrayList<double> results;
+ * double *params;
+ *
+ * mog.MakeModel(number_of_gaussians, dimension, params);
+ * mog.L2Error(data);
+ * mog.OutputResults(&results);
+ * @endcode
+ */
+class MoGL2E {
+ private:
+ // The parameters of the Mixture model
+ std::vector<arma::vec> mu_;
+ std::vector<arma::mat> sigma_;
+ arma::vec omega_;
+ size_t number_of_gaussians_;
+ size_t dimension_;
+
+ // The differential for the paramterization
+ // for optimization
+ arma::mat d_omega_;
+ std::vector<std::vector<arma::mat> > d_sigma_;
+
+ public:
+
+ MoGL2E() { }
+
+ ~MoGL2E() { }
+
+ void Init(size_t num_gauss, size_t dimension) {
+ // Destruct everything to initialize afresh
+ mu_.clear();
+ sigma_.clear();
+ d_sigma_.clear();
+
+ // Initialize the private variables
+ number_of_gaussians_ = num_gauss;
+ dimension_ = dimension;
+
+ // Resize the vector of vectors and matrices
+ mu_.resize(number_of_gaussians_);
+ sigma_.resize(number_of_gaussians_);
+ }
+
+ void Resize_d_sigma_() {
+ d_sigma_.resize(number_of_gaussians_);
+ for(size_t i = 0; i < number_of_gaussians_; i++)
+ d_sigma_[i].resize(dimension_ * (dimension_ + 1) / 2);
+ }
+
+ /**
+ * This function uses the parameters used for optimization
+ * and converts it into athe parameters of a Gaussian
+ * mixture model. This is to be used when you do not want
+ * the gradient values.
+ *
+ * Example use:
+ *
+ * @code
+ * MoGL2E mog;
+ * mog.MakeModel(number_of_gaussians, dimension,
+ * parameters_for_optimization);
+ * @endcode
+ */
+
+ void MakeModel(size_t num_mods, size_t dimension, const arma::vec& theta) {
+ arma::vec temp_mu(dimension);
+ arma::mat lower_triangle_matrix;
+ double sum, s_min = 0.01;
+
+ Init(num_mods, dimension);
+ lower_triangle_matrix.set_size(dimension, dimension);
+
+ // calculating the omega values
+ arma::vec temp_array = exp(theta);
+ temp_array[num_mods - 1] = 1;
+ sum = accu(temp_array);
+
+ temp_array /= sum;
+ set_omega(temp_array);
+
+ // calculating the mu values
+ for(size_t k = 0; k < num_mods; k++) {
+ for(size_t j = 0; j < dimension; j++) {
+ temp_mu[j] = theta[num_mods + (k * dimension) + j - 1];
+ }
+ set_mu(k, temp_mu);
+ }
+
+ // calculating the sigma values
+ // using a lower triangular matrix and its transpose
+ // to obtain a positive definite symmetric matrix
+ arma::mat sigma_temp(dimension, dimension);
+ for(size_t k = 0; k < num_mods; k++) {
+ lower_triangle_matrix.zeros();
+ for(size_t j = 0; j < dimension; j++) {
+ for(size_t i = 0; i < j; i++) {
+ lower_triangle_matrix(j, i) = theta[(num_mods - 1)
+ + (num_mods * dimension) + k * (dimension * (dimension + 1) / 2)
+ + (j * (j + 1) / 2) + i];
+ }
+ lower_triangle_matrix(j, j) = theta[(num_mods - 1)
+ + (num_mods * dimension) + k * (dimension * (dimension + 1) / 2)
+ + (j * (j + 1) / 2) + j] + s_min;
+ }
+ sigma_temp = lower_triangle_matrix * trans(lower_triangle_matrix);
+ set_sigma(k, sigma_temp);
+ }
+ }
+
+ /**
+ * This function uses the parameters used for optimization
+ * and converts it into the parameters of a Gaussian
+ * mixture model. This is to be used when you want
+ * the gradient values.
+ *
+ * Example use:
+ *
+ * @code
+ * MoGL2E mog;
+ * mog.MakeModelWithGradients(number_of_gaussians, dimension,
+ * parameters_for_optimization);
+ * @endcode
+ */
+ void MakeModelWithGradients(size_t num_mods,
+ size_t dimension,
+ const arma::vec& theta) {
+ arma::vec temp_mu(dimension);
+ arma::mat lower_triangle_matrix;
+ double sum, s_min = 0.01;
+
+ Init(num_mods, dimension);
+ lower_triangle_matrix.set_size(dimension, dimension);
+
+ // calculating the omega values
+ arma::vec temp_array = exp(theta);
+ temp_array[num_mods - 1] = 1;
+ sum = accu(temp_array);
+
+ temp_array /= sum;
+ set_omega(temp_array);
+
+ // calculating the d_omega values
+ arma::mat d_omega_temp(num_mods - 1, num_mods);
+ d_omega_temp.zeros();
+ for (size_t i = 0; i < num_mods - 1; i++) {
+ for (size_t j = 0; j < i; j++) {
+ d_omega_temp(i, j) = -(omega_[i] * omega_[j]);
+ d_omega_temp(j, i) = -(omega_[i] * omega_[j]);
+ }
+ d_omega_temp(i, i) = omega_[i] * (1 - omega_[i]);
+ }
+
+ for (size_t i = 0; i < num_mods - 1; i++)
+ d_omega_temp(i, num_mods - 1) = -(omega_[i] * omega_[num_mods - 1]);
+
+ set_d_omega(d_omega_temp);
+
+ // calculating the mu values
+ for (size_t k = 0; k < num_mods; k++) {
+ for (size_t j = 0; j < dimension; j++)
+ temp_mu[j] = theta[num_mods + (k * dimension) + j - 1];
+ set_mu(k, temp_mu);
+ }
+ // d_mu is not computed because it is implicitly known
+ // since no parameterization is applied on them
+
+ // using a lower triangular matrix and its transpose
+ // to obtain a positive definite symmetric matrix
+
+ // initializing the d_sigma values
+
+ arma::mat d_sigma_temp(dimension, dimension);
+ Resize_d_sigma_();
+
+ // calculating the sigma values
+ arma::mat sigma_temp(dimension, dimension);
+ for (size_t k = 0; k < num_mods; k++) {
+ lower_triangle_matrix.zeros();
+ for (size_t j = 0; j < dimension; j++) {
+ for (size_t i = 0; i < j; i++) {
+ lower_triangle_matrix(j, i) = theta[(num_mods - 1)
+ + (num_mods * dimension) + k * (dimension * (dimension + 1) / 2)
+ + (j * (j + 1) / 2) + i];
+ }
+ lower_triangle_matrix(j, j) = theta[(num_mods - 1)
+ + (num_mods * dimension) + k * (dimension * (dimension + 1) / 2)
+ + (j * (j + 1) / 2) + j] + s_min;
+ }
+ sigma_temp = lower_triangle_matrix * trans(lower_triangle_matrix);
+ set_sigma(k, sigma_temp);
+
+ // calculating the d_sigma values
+ for (size_t i = 0; i < dimension; i++) {
+ for (size_t in = 0; in < i + 1; in++) {
+ arma::mat d_sigma_d_r(dimension, dimension);
+ d_sigma_d_r.zeros();
+ d_sigma_d_r(i, in) = 1.0;
+
+ d_sigma_temp = d_sigma_d_r * trans(lower_triangle_matrix) +
+ lower_triangle_matrix * trans(d_sigma_d_r);
+ set_d_sigma(k, (i * (i + 1) / 2) + in, d_sigma_temp);
+ }
+ }
+ }
+ }
+
+ ////// THE GET FUNCTCLINS //////
+ std::vector<arma::vec>& mu() {
+ return mu_;
+ }
+
+ std::vector<arma::mat>& sigma() {
+ return sigma_;
+ }
+
+ arma::vec& omega() {
+ return omega_;
+ }
+
+ size_t number_of_gaussians() {
+ return number_of_gaussians_;
+ }
+
+ size_t dimension() {
+ return dimension_;
+ }
+
+ arma::vec& mu(size_t i) {
+ return mu_[i] ;
+ }
+
+ arma::mat& sigma(size_t i) {
+ return sigma_[i];
+ }
+
+ double omega(size_t i) {
+ return omega_[i];
+ }
+
+ arma::mat& d_omega() {
+ return d_omega_;
+ }
+
+ std::vector<std::vector<arma::mat> >& d_sigma() {
+ return d_sigma_;
+ }
+
+ std::vector<arma::mat>& d_sigma(size_t i) {
+ return d_sigma_[i];
+ }
+
+ ////// THE SET FUNCTCLINS //////
+
+ void set_mu(size_t i, const arma::vec& mu) {
+ assert(i < number_of_gaussians_);
+ assert(mu.n_elem == dimension_);
+
+ mu_[i] = mu;
+ }
+
+ void set_sigma(size_t i, const arma::mat& sigma) {
+ assert(i < number_of_gaussians_);
+ assert(sigma.n_rows == dimension_);
+ assert(sigma.n_cols == dimension_);
+
+ sigma_[i] = sigma;
+ }
+
+ void set_omega(const arma::vec& omega) {
+ assert(omega.n_elem == number_of_gaussians_);
+ omega_ = omega;
+ }
+
+ void set_d_omega(const arma::mat& d_omega) {
+ d_omega_ = d_omega;
+ }
+
+ void set_d_sigma(size_t i, size_t j, const arma::mat& d_sigma_i_j) {
+ d_sigma_[i][j] = d_sigma_i_j;
+ }
+
+ /**
+ * This function outputs the parameters of the model
+ * to an arraylist of doubles
+ *
+ * @code
+ * ArrayList<double> results;
+ * mog.OutputResults(&results);
+ * @endcode
+ */
+ void OutputResults(std::vector<double>& results) {
+
+ // Initialize the size of the output array
+ results.resize(number_of_gaussians_ * (1 + dimension_ * (1 + dimension_)));
+
+ // Copy values to the array from the private variables of the class
+ for (size_t i = 0; i < number_of_gaussians_; i++) {
+ results[i] = omega_[i];
+ for (size_t j = 0; j < dimension_; j++) {
+ results[number_of_gaussians_ + (i * dimension_) + j] = (mu_[i])[j];
+ for (size_t k = 0; k < dimension_; k++) {
+ results[number_of_gaussians_ * (1 + dimension_)
+ + (i * dimension_ * dimension_) + (j * dimension_)
+ + k] = (sigma_[i])(j, k);
+ }
+ }
+ }
+ }
+
+ /**
+ * This function prints the parameters of the model
+ *
+ * @code
+ * mog.Display();
+ * @endcode
+ */
+ void Display() {
+ // Output the model parameters as the omega, mu and sigma
+ Log::Info << " Omega : [ ";
+ for (size_t i = 0; i < number_of_gaussians_; i++) {
+ Log::Info << omega_[i];
+ }
+ Log::Info << "]" << std::endl;
+ Log::Info << " Mu : " << std::endl << "[";
+ for (size_t i = 0; i < number_of_gaussians_; i++) {
+ for (size_t j = 0; j < dimension_ ; j++) {
+ Log::Info << (mu_[i])[j];
+ }
+ Log::Info << ";";
+ if (i == (number_of_gaussians_ - 1)) {
+ Log::Info << "\b]" << std::endl;
+ }
+ }
+ Log::Info << "Sigma : ";
+ for (size_t i = 0; i < number_of_gaussians_; i++) {
+ Log::Info << std::endl << "[";
+ for (size_t j = 0; j < dimension_ ; j++) {
+ for (size_t k = 0; k < dimension_ ; k++) {
+ Log::Info << (sigma_[i])(j, k);
+ }
+ Log::Info << ";";
+ }
+ Log::Info << "\b]";
+ }
+ Log::Info << std::endl;
+ }
+
+ /**
+ * This function calculates the L2 error and
+ * the gradient of the error with respect to the
+ * parameters given the data and the parameterized
+ * mixture
+ *
+ * Example use:
+ *
+ * @code
+ * const Matrix data;
+ * MoGL2E mog;
+ * size_t num_gauss, dimension;
+ * double *params; // get the parameters
+ *
+ * mog.MakeModel(num_gauss, dimension, params);
+ * mog.L2Error(data);
+ * @endcode
+ */
+ long double L2Error(const arma::mat& data);
+ long double L2Error(const arma::mat& data, arma::vec& gradients);
+
+ /**
+ * Calculates the regularization value for a
+ * Gaussian mixture and its gradient with
+ * respect to the parameters
+ *
+ * Used by the 'L2Error' function to calculate
+ * the regularization part of the error
+ */
+ long double RegularizationTerm_();
+ long double RegularizationTerm_(arma::vec& g_reg);
+
+ /**
+ * Calculates the goodness-of-fit value for a
+ * Gaussian mixture and its gradient with
+ * respect to the parameters
+ *
+ * Used by the 'L2Error' function to calculate
+ * the goodness-of-fit part of the error
+ */
+ long double GoodnessOfFitTerm_(const arma::mat& data);
+ long double GoodnessOfFitTerm_(const arma::mat& data, arma::vec& g_fit);
+
+ /**
+ * This function computes multiple number of starting points
+ * required for the Nelder Mead method
+ *
+ * Example use:
+ * @code
+ * double **p;
+ * size_t n, num_gauss;
+ * const Matrix data;
+ *
+ * MoGL2E::MultiplePointsGeneratot(p, n, data, num_gauss);
+ * @endcode
+ */
+ static void MultiplePointsGenerator(arma::mat& points,
+ const arma::mat& d,
+ size_t number_of_components);
+
+ /**
+ * This function parameterizes the starting point obtained
+ * from the 'k_means" for optimization purposes using the
+ * Quasi Newton method
+ *
+ * Example use:
+ * @code
+ * double *p;
+ * size_t num_gauss;
+ * const Matrix data;
+ *
+ * MoGL2E::InitialPointGeneratot(p, data, num_gauss);
+ * @endcode
+ */
+ static void InitialPointGenerator(arma::vec& theta,
+ const arma::mat& data,
+ size_t k_comp);
+
+ /**
+ * This is the function which would be used for
+ * optimization. It creates its own object of
+ * class MoGL2E and returns the L2 error
+ * and the gradient which are computed by
+ * the functions of the class
+ *
+ */
+ static long double L2ErrorForOpt(const arma::vec& params,
+ const arma::mat& data) {
+ MoGL2E model;
+ size_t num_gauss = (params.n_elem + 1) * 2 /
+ ((data.n_rows + 1) * (data.n_rows + 2));
+
+ model.MakeModel(num_gauss, data.n_rows, params);
+
+ return model.L2Error(data);
+ }
+
+ static long double L2ErrorForOpt(const arma::vec& params,
+ const arma::mat& data,
+ arma::vec& gradient) {
+
+ MoGL2E model;
+ size_t num_gauss = (params.n_elem + 1) * 2 /
+ ((data.n_rows + 1) * (data.n_rows + 2));
+
+ model.MakeModelWithGradients(num_gauss, data.n_rows, params);
+
+ return model.L2Error(data, gradient);
+ }
+};
+
+}; // namespace gmm
+}; // namespace mlpack
+
+#endif // __MLPACK_METHODS_MOG_MOG_L2E_HPP
Deleted: mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cc
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cc 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cc 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,122 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_l2e_main.cc
- *
- * This program test drives the L2 estimation
- * of a Gaussian Mixture model.
- *
- * PARAMETERS TO BE INPUT:
- *
- * --data
- * This is the file that contains the data on which
- * the model is to be fit
- *
- * --mog_l2e/K
- * This is the number of gaussians we want to fit
- * on the data, defaults to '1'
- *
- * --output
- * This file will contain the parameters estimated,
- * defaults to 'output.csv'
- *
- */
-#include "mog_l2e.h"
-#include "optimizers.h"
-
-PROGRAM_INFO("Mixture of Gaussians",
- "This program takes a parametric estimate of a Gaussian mixture model (GMM)"
- " using the L2 loss function.", "mog_l2e");
-
-PARAM_STRING_REQ("data", "A file containing the data on which the model has to "
- "be fit.", "mog_l2e");
-PARAM_STRING("output", "The file into which the output is to be written into",
- "mog_l2e", "output.csv");
-
-using namespace mlpack;
-using namespace mlpack::gmm;
-
-int main(int argc, char* argv[]) {
- CLI::ParseCommandLine(argc, argv);
-
- ////// READING PARAMETERS AND LOADING DATA //////
- arma::mat data_points;
- data_points.load(CLI::GetParam<std::string>("mog_l2e/data").c_str(),
- arma::auto_detect, false, true);
-
- ////// MIXTURE OF GAUSSIANS USING L2 ESTIMATCLIN //////
- size_t number_of_gaussians = CLI::GetParam<int>("mog_l2e/k");
- CLI::GetParam<int>("mog_l2e/d") = data_points.n_rows;
- size_t dimension = CLI::GetParam<int>("mog_l2e/d");
-
- ////// RUNNING AN OPTIMIZER TO MINIMIZE THE L2 ERROR //////
- const char *opt_method = CLI::GetParam<std::string>("opt/method").c_str();
- size_t param_dim = (number_of_gaussians * (dimension + 1) * (dimension + 2)
- / 2 - 1);
- CLI::GetParam<int>("opt/param_space_dim") = param_dim;
-
- size_t optim_flag = (strcmp(opt_method, "NelderMead") == 0 ? 1 : 0);
- MoGL2E mog;
-
- if (optim_flag == 1) {
- ////// OPTIMIZER USING NELDER MEAD METHOD //////
- NelderMead opt;
-
- ////// Initializing the optimizer //////
- CLI::StartTimer("opt/init_opt");
- opt.Init(MoGL2E::L2ErrorForOpt, data_points);
- CLI::StopTimer("opt/init_opt");
-
- ////// Getting starting points for the optimization //////
- arma::mat pts(param_dim, param_dim + 1);
-
- CLI::StartTimer("opt/get_init_pts");
- MoGL2E::MultiplePointsGenerator(pts, data_points, number_of_gaussians);
- CLI::StopTimer("opt/get_init_pts");
-
- ////// The optimization //////
- CLI::StartTimer("opt/optimizing");
- opt.Eval(pts);
- CLI::StopTimer("opt/optimizing");
-
- ////// Making model with the optimal parameters //////
- // This is a stupid way to do it and putting the 0s there ensures it will
- // fail. Do it better!
- mog.MakeModel(0, 0, pts.col(0));
- } else {
- ////// OPTIMIZER USING QUASI NEWTON METHOD //////
- QuasiNewton opt;
-
- ////// Initializing the optimizer //////
- CLI::StartTimer("opt/init_opt");
- opt.Init(MoGL2E::L2ErrorForOpt, data_points);
- CLI::StopTimer("opt/init_opt");
-
- ////// Getting starting point for the optimization //////
- arma::vec pt(param_dim);
-
- CLI::StartTimer("opt/get_init_pt");
- MoGL2E::InitialPointGenerator(pt, data_points, number_of_gaussians);
- CLI::StopTimer("opt/get_init_pt");
-
- ////// The optimization //////
- CLI::StartTimer("opt/optimizing");
- opt.Eval(pt);
- CLI::StopTimer("opt/optimizing");
-
- ////// Making model with optimal parameters //////
- // This is a stupid way to do it and putting the 0s there ensures it will
- // fail. Do it better!
- mog.MakeModel(0, 0, pt);
-
- }
-
- long double error = mog.L2Error(data_points);
- Log::Info << "Minimum L2 error achieved: " << error << "." << std::endl;
- mog.Display();
-
- std::vector<double> results;
- mog.OutputResults(results);
-
- ////// OUTPUT RESULTS //////
- // We need a better way to do this (like XML). For now we just do nothing.
-}
Copied: mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cc)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/mog_l2e_main.cpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,122 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_l2e_main.cpp
+ *
+ * This program test drives the L2 estimation
+ * of a Gaussian Mixture model.
+ *
+ * PARAMETERS TO BE INPUT:
+ *
+ * --data
+ * This is the file that contains the data on which
+ * the model is to be fit
+ *
+ * --mog_l2e/K
+ * This is the number of gaussians we want to fit
+ * on the data, defaults to '1'
+ *
+ * --output
+ * This file will contain the parameters estimated,
+ * defaults to 'output.csv'
+ *
+ */
+#include "mog_l2e.hpp"
+#include "optimizers.hpp"
+
+PROGRAM_INFO("Mixture of Gaussians",
+ "This program takes a parametric estimate of a Gaussian mixture model (GMM)"
+ " using the L2 loss function.", "mog_l2e");
+
+PARAM_STRING_REQ("data", "A file containing the data on which the model has to "
+ "be fit.", "mog_l2e");
+PARAM_STRING("output", "The file into which the output is to be written into",
+ "mog_l2e", "output.csv");
+
+using namespace mlpack;
+using namespace mlpack::gmm;
+
+int main(int argc, char* argv[]) {
+ CLI::ParseCommandLine(argc, argv);
+
+ ////// READING PARAMETERS AND LOADING DATA //////
+ arma::mat data_points;
+ data_points.load(CLI::GetParam<std::string>("mog_l2e/data").c_str(),
+ arma::auto_detect, false, true);
+
+ ////// MIXTURE OF GAUSSIANS USING L2 ESTIMATCLIN //////
+ size_t number_of_gaussians = CLI::GetParam<int>("mog_l2e/k");
+ CLI::GetParam<int>("mog_l2e/d") = data_points.n_rows;
+ size_t dimension = CLI::GetParam<int>("mog_l2e/d");
+
+ ////// RUNNING AN OPTIMIZER TO MINIMIZE THE L2 ERROR //////
+ const char *opt_method = CLI::GetParam<std::string>("opt/method").c_str();
+ size_t param_dim = (number_of_gaussians * (dimension + 1) * (dimension + 2)
+ / 2 - 1);
+ CLI::GetParam<int>("opt/param_space_dim") = param_dim;
+
+ size_t optim_flag = (strcmp(opt_method, "NelderMead") == 0 ? 1 : 0);
+ MoGL2E mog;
+
+ if (optim_flag == 1) {
+ ////// OPTIMIZER USING NELDER MEAD METHOD //////
+ NelderMead opt;
+
+ ////// Initializing the optimizer //////
+ CLI::StartTimer("opt/init_opt");
+ opt.Init(MoGL2E::L2ErrorForOpt, data_points);
+ CLI::StopTimer("opt/init_opt");
+
+ ////// Getting starting points for the optimization //////
+ arma::mat pts(param_dim, param_dim + 1);
+
+ CLI::StartTimer("opt/get_init_pts");
+ MoGL2E::MultiplePointsGenerator(pts, data_points, number_of_gaussians);
+ CLI::StopTimer("opt/get_init_pts");
+
+ ////// The optimization //////
+ CLI::StartTimer("opt/optimizing");
+ opt.Eval(pts);
+ CLI::StopTimer("opt/optimizing");
+
+ ////// Making model with the optimal parameters //////
+ // This is a stupid way to do it and putting the 0s there ensures it will
+ // fail. Do it better!
+ mog.MakeModel(0, 0, pts.col(0));
+ } else {
+ ////// OPTIMIZER USING QUASI NEWTON METHOD //////
+ QuasiNewton opt;
+
+ ////// Initializing the optimizer //////
+ CLI::StartTimer("opt/init_opt");
+ opt.Init(MoGL2E::L2ErrorForOpt, data_points);
+ CLI::StopTimer("opt/init_opt");
+
+ ////// Getting starting point for the optimization //////
+ arma::vec pt(param_dim);
+
+ CLI::StartTimer("opt/get_init_pt");
+ MoGL2E::InitialPointGenerator(pt, data_points, number_of_gaussians);
+ CLI::StopTimer("opt/get_init_pt");
+
+ ////// The optimization //////
+ CLI::StartTimer("opt/optimizing");
+ opt.Eval(pt);
+ CLI::StopTimer("opt/optimizing");
+
+ ////// Making model with optimal parameters //////
+ // This is a stupid way to do it and putting the 0s there ensures it will
+ // fail. Do it better!
+ mog.MakeModel(0, 0, pt);
+
+ }
+
+ long double error = mog.L2Error(data_points);
+ Log::Info << "Minimum L2 error achieved: " << error << "." << std::endl;
+ mog.Display();
+
+ std::vector<double> results;
+ mog.OutputResults(results);
+
+ ////// OUTPUT RESULTS //////
+ // We need a better way to do this (like XML). For now we just do nothing.
+}
Deleted: mlpack/trunk/src/mlpack/methods/mog/optimizers.cc
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/optimizers.cc 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/optimizers.cc 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,351 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file optimizers.cc
- *
- * Implementation of the optimizers
- */
-#include <mlpack/core.h>
-#include "optimizers.h"
-
-using namespace mlpack;
-
-void NelderMead::Eval(arma::mat& pts) {
- size_t dim = dimension();
- size_t num_func_eval;
- size_t i, j, ihi, ilo, inhi;
- size_t mpts = dim + 1;
- double sum, swap;
- arma::vec psum;
- long double swap_y, rtol, ytry, TINY = 1.0e-10;
- arma::vec y;
- arma::vec param_passed;
-
- // Default value is 1e-5.
- long double tol = CLI::GetParam<double>("opt/tolerance");
-
- // Default value is 50000.
- size_t NMAX = CLI::GetParam<int>("opt/MAX_FUNC_EVAL");
-
- param_passed.set_size(dim);
- psum.set_size(dim);
- num_func_eval = 0;
- y.set_size(mpts);
- for(i = 0; i < mpts; i++) {
- param_passed = pts.row(i);
- y[i] = (*func_ptr_)(param_passed, data_);
- }
-
- for(;;) {
- ilo = 0;
- ihi = y[0] > y[1] ? (inhi = 1, 0) : (inhi = 0, 1);
- for (i = 0; i < mpts; i++) {
- if(y[i] <= y[ilo])
- ilo = i;
- if(y[i] > y[ihi]) {
- inhi = ihi;
- ihi = i;
- }
- else if((y[i] > y[inhi]) && (i != ihi))
- inhi = i;
- }
-
- rtol = 2.0 * fabs(y[ihi] - y[ilo]) / (fabs(y[ihi]) + fabs(y[ilo]) + TINY);
- if (rtol < tol) {
- swap_y = y[0];
- y[0] = y[ilo];
- y[ilo] = swap_y;
- for (i = 0; i < dim; i++) {
- swap = pts(0, i);
- pts(0, i) = pts(ilo, i);
- pts(ilo, i) = swap;
- }
-
- break;
- }
-
- if (num_func_eval > NMAX) {
- Log::Warn << "Nelder-Mead: Maximum number of function evaluations "
- "exceeded." << std::endl;
- break;
- }
-
- num_func_eval += 2;
-
- // Beginning a new iteration.
- // Extrapolating by a factor of -1.0 through the face of the simplex
- // across from the high point, i.e, reflect the simplex from the high point
- for (j = 0 ; j < dim; j++) {
- sum = 0.0;
- for(i = 0; i < mpts; i++)
- if (i != ihi)
- sum += pts(i, j);
-
- psum[j] = sum / dim;
- }
-
- ytry = ModSimplex_(pts, y, psum, ihi, -1.0);
-
- if(ytry <= y[ilo]) {
- // result better than best point
- // so additional extrapolation by a factor of 2
- ytry = ModSimplex_(pts, y, psum, ihi, 2.0);
- } else if(ytry >= y[ihi]) {
- // result worse than the worst point
- // so there is a lower intermediate point,
- // i.e., do a one dimensional contraction
- ytry = ModSimplex_(pts, y, psum, ihi, 0.5);
- if(ytry > y[ihi]) {
- // Can't get rid of the high point,
- // try to contract around the best point
- for (i = 0; i < mpts; i++) {
- if (i != ilo) {
- for (j = 0; j < dim; j++)
- pts(i, j) = psum[j] = 0.5 * (pts(i, j) + pts(ilo, j));
-
- param_passed = psum;
- y[i] = (*func_ptr_)(param_passed, data());
- }
- }
- num_func_eval += dim;
-
- for (j = 0 ; j < dim ; j++) {
- sum = 0.0;
- for (i = 0 ; i < mpts ; i++)
- if (i != ihi)
- sum += pts(i, j);
-
- psum[j] = sum / dim;
- }
- }
- }
- else
- --num_func_eval;
- }
-
- CLI::GetParam<int>("opt/func_evals") = num_func_eval;
-}
-
-long double NelderMead::ModSimplex_(arma::mat& pts, arma::vec& y,
- arma::vec& psum, size_t ihi,
- float fac) {
-
- size_t j, dim = dimension();
- long double ytry;
- arma::vec ptry(dim);
- arma::vec param_passed(dim);
-
- for (j = 0; j < dim; j++)
- ptry[j] = psum[j] * (1 - fac) + pts(ihi, j) * fac;
-
- param_passed = ptry;
- ytry = (*func_ptr_)(param_passed, data());
- if (ytry < y[ihi]) {
- y[ihi] = ytry;
- for (j = 0; j < dim; j++)
- pts(ihi, j) = ptry[j];
- }
-
- return ytry;
-}
-
-void QuasiNewton::Eval(arma::vec& pt) {
- size_t n = dimension(), iters;
- size_t i, its;
- // Default value is 200.
- size_t MAXIMUM_ITERATCLINS = CLI::GetParam<int>("opt/MAX_ITERS");
-
- long double temp_1, temp_2, temp_3, temp_4, f_previous, f_min,
- maximum_step_length, sum = 0.0, sumdg, sumxi, temp, test;
- arma::vec dgrad, grad, hdgrad, xi;
- arma::vec pold, pnew;
- arma::mat hessian;
-
- // Default value is 3.0e-8.
- double EPSILON = CLI::GetParam<double>("opt/EPSILON");
-
- // Default value is 1.0e-5.
- double TOLERANCE = CLI::GetParam<double>("optTOLERANCE");
-
- // Default value is 100.
- double MAX_STEP_SIZE = CLI::GetParam<double>("opt/MAX_STEP_SIZE");
-
- // Default value is 1.0e-7.
- double g_tol = CLI::GetParam<double>("opt/gtol");
-
- dgrad.set_size(n);
- grad.set_size(n);
- hdgrad.set_size(n);
- hessian.set_size(n, n);
- pnew.set_size(n);
- xi.set_size(n);
- pold = pt;
- f_previous = (*func_ptr_)(pold, data(), grad);
- arma::vec tmp;
- tmp.ones(n);
- hessian.diag() = tmp;
-
- xi = -1 * grad;
- sum = dot(pold, pold);
-
- double fmax;
- if(sqrt(sum) > (float) n)
- fmax = sqrt(sum);
- else
- fmax = (float) n;
-
- maximum_step_length = MAX_STEP_SIZE * fmax;
-
- for (its = 0; its < MAXIMUM_ITERATCLINS; its++) {
- dgrad = grad;
- LineSearch_(pold, f_previous, grad, xi,
- pnew, f_min, maximum_step_length);
- f_previous = f_min;
- xi = pold - pnew;
- pold = pnew;
- pt = pold;
-
- test = 0.0;
- for (i = 0; i < n; i++) {
- if (fabs(pold[i]) > 1.0)
- fmax = fabs(pold[i]);
- else
- fmax = 1.0;
-
- temp = fabs(xi[i]) / fmax;
- if (temp > test)
- test = temp;
- }
-
- if (test < TOLERANCE) {
- iters = its;
- CLI::GetParam<int>("opt/iters") = iters;
- return;
- }
-
- test = 0.0;
- if (f_min > 1.0)
- temp_1 = f_min;
- else
- temp_1 = 1.0;
-
- for (i = 0; i < n; i++) {
- if (fabs(pold[i]) > 1.0)
- fmax = pold[i];
- else
- fmax = 1.0;
-
- temp = fabs(grad[i]) * fmax / temp_1;
- if (temp > test)
- test = temp;
- }
-
- if (test < g_tol) {
- iters = its;
- CLI::GetParam<int>("opt/iters") = iters;
- return;
- }
-
- dgrad -= grad;
- dgrad *= -1.0;
- hdgrad = hessian * dgrad;
-
- temp_2 = dot(dgrad, xi);
- temp_4 = dot(dgrad, hdgrad);
- sumdg = dot(dgrad, dgrad);
- sumxi = dot(xi, xi);
-
- if (temp_2 > sqrt(EPSILON * sumdg * sumxi)) {
- temp_2 = 1.0 / temp_2;
- temp_3 = 1.0 / temp_4;
-
- dgrad = temp_2 * xi;
- dgrad -= temp_3 * hdgrad;
-
- hessian += temp_2 * (xi * trans(xi));
- hessian -= temp_3 * (hdgrad * trans(hdgrad));
- hessian += temp_4 * (dgrad * trans(dgrad));
- }
-
- xi = -hessian * grad;
- }
-
- Log::Warn << "Too many iterations in Quasi-Newton: giving up." << std::endl;
-}
-
-void QuasiNewton::LineSearch_(arma::vec& pold, long double fold,
- arma::vec& grad, arma::vec& xi,
- arma::vec& pnew, long double& f_min,
- long double maximum_step_length) {
-
- size_t i, n = dimension();
- long double a, step_length, previous_step_length = 0.0,
- minimum_step_length, b, disc, previous_f_value = 0.0,
- rhs1, rhs2, slope, sum, temp, test, temp_step_length,
- MIN_DECREASE = 1.0e-4, TOLERANCE = 1.0e-7;
-
- sum = sqrt(dot(xi, xi));
- if(sum > maximum_step_length)
- xi *= (maximum_step_length / sum);
-
- slope = dot(grad, xi);
- if (slope >= 0.0)
- return;
-
- test = 0.0;
- for (i = 0; i < n; i++) {
- double fmax;
- fmax = (fabs(pold[i]) > 1.0 ? fabs(pold[i]) : 1.0);
- temp = fabs(xi[i]) / fmax;
- if (temp > test)
- test = temp;
- }
-
- minimum_step_length = TOLERANCE / test;
- step_length = 1.0;
-
- for (;;) {
- pnew = pold + step_length * xi;
-
- f_min = (*func_ptr_)(pnew, data(), grad);
-
- if (step_length < minimum_step_length) {
- pnew = pold;
- return;
- } else if(f_min <= fold + MIN_DECREASE * step_length * slope) {
- return;
- } else {
- if (step_length == 1.0)
- temp_step_length = -slope / (2.0 * (f_min - fold - slope));
- else {
- rhs1 = f_min - fold - step_length * slope;
- rhs2 = previous_f_value - fold - previous_step_length * slope;
- a = (rhs1 / (step_length * step_length)
- - rhs2 / (previous_step_length * previous_step_length))
- / (step_length - previous_step_length);
- b = (-previous_step_length * rhs1 / (step_length * step_length) +
- step_length * rhs2 / (previous_step_length * previous_step_length))
- / (step_length - previous_step_length);
-
- if (a == 0.0)
- temp_step_length = -slope / (2.0*b);
- else {
- disc = b * b - 3.0 * a * slope;
- if(disc < 0.0)
- temp_step_length = 0.5 * step_length;
- else if (b <= 0.0)
- temp_step_length = (-b + sqrt(disc)) / (3.0 * a);
- else
- temp_step_length = -slope / (b + sqrt(disc));
- }
-
- if(temp_step_length > 0.5 * step_length)
- temp_step_length = 0.5 * step_length;
- }
- }
-
- previous_step_length = step_length;
- previous_f_value = f_min;
- step_length = (temp_step_length > 0.1 * step_length) ? temp_step_length
- : 0.1 * step_length;
- }
-}
Copied: mlpack/trunk/src/mlpack/methods/mog/optimizers.cpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/optimizers.cc)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/optimizers.cpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/optimizers.cpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,352 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file optimizers.cpp
+ *
+ * Implementation of the optimizers
+ */
+#include <mlpack/core.h>
+#include "optimizers.hpp"
+
+using namespace mlpack;
+using namespace gmm;
+
+void NelderMead::Eval(arma::mat& pts) {
+ size_t dim = dimension();
+ size_t num_func_eval;
+ size_t i, j, ihi, ilo, inhi;
+ size_t mpts = dim + 1;
+ double sum, swap;
+ arma::vec psum;
+ long double swap_y, rtol, ytry, TINY = 1.0e-10;
+ arma::vec y;
+ arma::vec param_passed;
+
+ // Default value is 1e-5.
+ long double tol = CLI::GetParam<double>("opt/tolerance");
+
+ // Default value is 50000.
+ size_t NMAX = CLI::GetParam<int>("opt/MAX_FUNC_EVAL");
+
+ param_passed.set_size(dim);
+ psum.set_size(dim);
+ num_func_eval = 0;
+ y.set_size(mpts);
+ for(i = 0; i < mpts; i++) {
+ param_passed = pts.row(i);
+ y[i] = (*func_ptr_)(param_passed, data_);
+ }
+
+ for(;;) {
+ ilo = 0;
+ ihi = y[0] > y[1] ? (inhi = 1, 0) : (inhi = 0, 1);
+ for (i = 0; i < mpts; i++) {
+ if(y[i] <= y[ilo])
+ ilo = i;
+ if(y[i] > y[ihi]) {
+ inhi = ihi;
+ ihi = i;
+ }
+ else if((y[i] > y[inhi]) && (i != ihi))
+ inhi = i;
+ }
+
+ rtol = 2.0 * fabs(y[ihi] - y[ilo]) / (fabs(y[ihi]) + fabs(y[ilo]) + TINY);
+ if (rtol < tol) {
+ swap_y = y[0];
+ y[0] = y[ilo];
+ y[ilo] = swap_y;
+ for (i = 0; i < dim; i++) {
+ swap = pts(0, i);
+ pts(0, i) = pts(ilo, i);
+ pts(ilo, i) = swap;
+ }
+
+ break;
+ }
+
+ if (num_func_eval > NMAX) {
+ Log::Warn << "Nelder-Mead: Maximum number of function evaluations "
+ "exceeded." << std::endl;
+ break;
+ }
+
+ num_func_eval += 2;
+
+ // Beginning a new iteration.
+ // Extrapolating by a factor of -1.0 through the face of the simplex
+ // across from the high point, i.e, reflect the simplex from the high point
+ for (j = 0 ; j < dim; j++) {
+ sum = 0.0;
+ for(i = 0; i < mpts; i++)
+ if (i != ihi)
+ sum += pts(i, j);
+
+ psum[j] = sum / dim;
+ }
+
+ ytry = ModSimplex_(pts, y, psum, ihi, -1.0);
+
+ if(ytry <= y[ilo]) {
+ // result better than best point
+ // so additional extrapolation by a factor of 2
+ ytry = ModSimplex_(pts, y, psum, ihi, 2.0);
+ } else if(ytry >= y[ihi]) {
+ // result worse than the worst point
+ // so there is a lower intermediate point,
+ // i.e., do a one dimensional contraction
+ ytry = ModSimplex_(pts, y, psum, ihi, 0.5);
+ if(ytry > y[ihi]) {
+ // Can't get rid of the high point,
+ // try to contract around the best point
+ for (i = 0; i < mpts; i++) {
+ if (i != ilo) {
+ for (j = 0; j < dim; j++)
+ pts(i, j) = psum[j] = 0.5 * (pts(i, j) + pts(ilo, j));
+
+ param_passed = psum;
+ y[i] = (*func_ptr_)(param_passed, data());
+ }
+ }
+ num_func_eval += dim;
+
+ for (j = 0 ; j < dim ; j++) {
+ sum = 0.0;
+ for (i = 0 ; i < mpts ; i++)
+ if (i != ihi)
+ sum += pts(i, j);
+
+ psum[j] = sum / dim;
+ }
+ }
+ }
+ else
+ --num_func_eval;
+ }
+
+ CLI::GetParam<int>("opt/func_evals") = num_func_eval;
+}
+
+long double NelderMead::ModSimplex_(arma::mat& pts, arma::vec& y,
+ arma::vec& psum, size_t ihi,
+ float fac) {
+
+ size_t j, dim = dimension();
+ long double ytry;
+ arma::vec ptry(dim);
+ arma::vec param_passed(dim);
+
+ for (j = 0; j < dim; j++)
+ ptry[j] = psum[j] * (1 - fac) + pts(ihi, j) * fac;
+
+ param_passed = ptry;
+ ytry = (*func_ptr_)(param_passed, data());
+ if (ytry < y[ihi]) {
+ y[ihi] = ytry;
+ for (j = 0; j < dim; j++)
+ pts(ihi, j) = ptry[j];
+ }
+
+ return ytry;
+}
+
+void QuasiNewton::Eval(arma::vec& pt) {
+ size_t n = dimension(), iters;
+ size_t i, its;
+ // Default value is 200.
+ size_t MAXIMUM_ITERATCLINS = CLI::GetParam<int>("opt/MAX_ITERS");
+
+ long double temp_1, temp_2, temp_3, temp_4, f_previous, f_min,
+ maximum_step_length, sum = 0.0, sumdg, sumxi, temp, test;
+ arma::vec dgrad, grad, hdgrad, xi;
+ arma::vec pold, pnew;
+ arma::mat hessian;
+
+ // Default value is 3.0e-8.
+ double EPSILON = CLI::GetParam<double>("opt/EPSILON");
+
+ // Default value is 1.0e-5.
+ double TOLERANCE = CLI::GetParam<double>("optTOLERANCE");
+
+ // Default value is 100.
+ double MAX_STEP_SIZE = CLI::GetParam<double>("opt/MAX_STEP_SIZE");
+
+ // Default value is 1.0e-7.
+ double g_tol = CLI::GetParam<double>("opt/gtol");
+
+ dgrad.set_size(n);
+ grad.set_size(n);
+ hdgrad.set_size(n);
+ hessian.set_size(n, n);
+ pnew.set_size(n);
+ xi.set_size(n);
+ pold = pt;
+ f_previous = (*func_ptr_)(pold, data(), grad);
+ arma::vec tmp;
+ tmp.ones(n);
+ hessian.diag() = tmp;
+
+ xi = -1 * grad;
+ sum = dot(pold, pold);
+
+ double fmax;
+ if(sqrt(sum) > (float) n)
+ fmax = sqrt(sum);
+ else
+ fmax = (float) n;
+
+ maximum_step_length = MAX_STEP_SIZE * fmax;
+
+ for (its = 0; its < MAXIMUM_ITERATCLINS; its++) {
+ dgrad = grad;
+ LineSearch_(pold, f_previous, grad, xi,
+ pnew, f_min, maximum_step_length);
+ f_previous = f_min;
+ xi = pold - pnew;
+ pold = pnew;
+ pt = pold;
+
+ test = 0.0;
+ for (i = 0; i < n; i++) {
+ if (fabs(pold[i]) > 1.0)
+ fmax = fabs(pold[i]);
+ else
+ fmax = 1.0;
+
+ temp = fabs(xi[i]) / fmax;
+ if (temp > test)
+ test = temp;
+ }
+
+ if (test < TOLERANCE) {
+ iters = its;
+ CLI::GetParam<int>("opt/iters") = iters;
+ return;
+ }
+
+ test = 0.0;
+ if (f_min > 1.0)
+ temp_1 = f_min;
+ else
+ temp_1 = 1.0;
+
+ for (i = 0; i < n; i++) {
+ if (fabs(pold[i]) > 1.0)
+ fmax = pold[i];
+ else
+ fmax = 1.0;
+
+ temp = fabs(grad[i]) * fmax / temp_1;
+ if (temp > test)
+ test = temp;
+ }
+
+ if (test < g_tol) {
+ iters = its;
+ CLI::GetParam<int>("opt/iters") = iters;
+ return;
+ }
+
+ dgrad -= grad;
+ dgrad *= -1.0;
+ hdgrad = hessian * dgrad;
+
+ temp_2 = dot(dgrad, xi);
+ temp_4 = dot(dgrad, hdgrad);
+ sumdg = dot(dgrad, dgrad);
+ sumxi = dot(xi, xi);
+
+ if (temp_2 > sqrt(EPSILON * sumdg * sumxi)) {
+ temp_2 = 1.0 / temp_2;
+ temp_3 = 1.0 / temp_4;
+
+ dgrad = temp_2 * xi;
+ dgrad -= temp_3 * hdgrad;
+
+ hessian += temp_2 * (xi * trans(xi));
+ hessian -= temp_3 * (hdgrad * trans(hdgrad));
+ hessian += temp_4 * (dgrad * trans(dgrad));
+ }
+
+ xi = -hessian * grad;
+ }
+
+ Log::Warn << "Too many iterations in Quasi-Newton: giving up." << std::endl;
+}
+
+void QuasiNewton::LineSearch_(arma::vec& pold, long double fold,
+ arma::vec& grad, arma::vec& xi,
+ arma::vec& pnew, long double& f_min,
+ long double maximum_step_length) {
+
+ size_t i, n = dimension();
+ long double a, step_length, previous_step_length = 0.0,
+ minimum_step_length, b, disc, previous_f_value = 0.0,
+ rhs1, rhs2, slope, sum, temp, test, temp_step_length,
+ MIN_DECREASE = 1.0e-4, TOLERANCE = 1.0e-7;
+
+ sum = sqrt(dot(xi, xi));
+ if(sum > maximum_step_length)
+ xi *= (maximum_step_length / sum);
+
+ slope = dot(grad, xi);
+ if (slope >= 0.0)
+ return;
+
+ test = 0.0;
+ for (i = 0; i < n; i++) {
+ double fmax;
+ fmax = (fabs(pold[i]) > 1.0 ? fabs(pold[i]) : 1.0);
+ temp = fabs(xi[i]) / fmax;
+ if (temp > test)
+ test = temp;
+ }
+
+ minimum_step_length = TOLERANCE / test;
+ step_length = 1.0;
+
+ for (;;) {
+ pnew = pold + step_length * xi;
+
+ f_min = (*func_ptr_)(pnew, data(), grad);
+
+ if (step_length < minimum_step_length) {
+ pnew = pold;
+ return;
+ } else if(f_min <= fold + MIN_DECREASE * step_length * slope) {
+ return;
+ } else {
+ if (step_length == 1.0)
+ temp_step_length = -slope / (2.0 * (f_min - fold - slope));
+ else {
+ rhs1 = f_min - fold - step_length * slope;
+ rhs2 = previous_f_value - fold - previous_step_length * slope;
+ a = (rhs1 / (step_length * step_length)
+ - rhs2 / (previous_step_length * previous_step_length))
+ / (step_length - previous_step_length);
+ b = (-previous_step_length * rhs1 / (step_length * step_length) +
+ step_length * rhs2 / (previous_step_length * previous_step_length))
+ / (step_length - previous_step_length);
+
+ if (a == 0.0)
+ temp_step_length = -slope / (2.0*b);
+ else {
+ disc = b * b - 3.0 * a * slope;
+ if(disc < 0.0)
+ temp_step_length = 0.5 * step_length;
+ else if (b <= 0.0)
+ temp_step_length = (-b + sqrt(disc)) / (3.0 * a);
+ else
+ temp_step_length = -slope / (b + sqrt(disc));
+ }
+
+ if(temp_step_length > 0.5 * step_length)
+ temp_step_length = 0.5 * step_length;
+ }
+ }
+
+ previous_step_length = step_length;
+ previous_f_value = f_min;
+ step_length = (temp_step_length > 0.1 * step_length) ? temp_step_length
+ : 0.1 * step_length;
+ }
+}
Deleted: mlpack/trunk/src/mlpack/methods/mog/optimizers.h
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/optimizers.h 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/optimizers.h 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,147 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file optimizers.h
- *
- * Declares classes for two types of optimizer
- *
- */
-#ifndef __MLPACK_METHODS_MOG_OPTIMIZERS_H
-#define __MLPACK_METHODS_MOG_OPTIMIZERS_H
-
-#include <mlpack/core.h>
-
-PARAM_STRING("method", "The method used to optimize", "opt", "");
-
-PARAM_INT_REQ("param_space_dim", "The dimension of the parameter space.", "opt");
-PARAM_INT("MAX_FUNC_EVAL", "The maximum number of function evaluations\
- allowed to the NelderMead optimizer (defaults to 50000)", "opt", 50000);
-
-PARAM_INT("func_evals", "The number of function evaluations taken by the algorithm", "opt", 0);
-PARAM_INT("MAX_ITERS", "The maximum number of iterations allowed to the function", "opt", 200);
-PARAM_INT("iters", "The number of iterations the algorithm actually went through", "opt", 0);
-
-PARAM_DOUBLE("EPSILON", "Value of epsilon.", "opt", 3.0e-8);
-PARAM_DOUBLE("TOLERANCE", "Tolerance for the minimum movement for the parameter value.", "opt", 1.0e-5);
-PARAM_DOUBLE("gtol", "Tolerance value for the gradient of the function", "opt", 1.0e-7);
-PARAM_DOUBLE("MAX_STEP_SIZE", "The maximum step size in the direction of the gradient.", "opt", 100.0);
-PARAM_DOUBLE("tolerance", "Undocumented parameter", "opt", 1.0e-5);
-
-PARAM_MODULE("opt", "This file contains two optimizers.");
-
-/**
- * An optimizer using the Nelder Mead method,
- * also known as the polytope or the simplex
- * method.
- *
- * It does multivariate minimization of an
- * objective function. If it is optimizing in
- * 'd' dimensions, it would require 'd+1'
- * starting points.
- *
- * Example use:
- *
- * @code
- * double init_pts[d+1][d];
- * size_t number_of_function_evaluations;
- * struct datanode *opt_module = fx_submodule(NULL,"NelderMead","opt_module");
- * Matrix data;
- * size_t dim_param_space;
- *
- * ...
- * NelderMead opt;
- * opt.Init(obj_function, data, dim_param_space, opt_module);
- * ...
- * opt.Eval(init_pts);
- * // init_pts[0] contaings the optimal point found
- * @endcode
- *
- */
-class NelderMead {
- private:
- size_t dimension_;
- arma::mat data_;
- long double (*func_ptr_)(const arma::vec&, const arma::mat&);
-
- public:
- NelderMead() { }
- ~NelderMead() { }
-
- void Init(long double (*fun)(const arma::vec&, const arma::mat&),
- arma::mat& data) {
- data_ = data;
- func_ptr_ = fun;
- dimension_ = mlpack::CLI::GetParam<int>("opt/param_space_dim");
- }
-
- const arma::mat& data() {
- return data_;
- }
-
- size_t dimension() {
- return dimension_;
- }
-
- void Eval(arma::mat& pts);
- long double ModSimplex_(arma::mat& pts, arma::vec& y,
- arma::vec& psum, size_t ihi, float fac);
-};
-
-/**
- * An optimizer using the Quasi Newton method,
- * also known as the variable metrics
- * method.
- *
- * It does multivariate minimization of an
- * objective function using only the function
- * value and the gradients.
- *
- * Example use:
- *
- * @code
- * double init_pt[d];
- * size_t number_of_iters;
- * struct datanode *opt_module = fx_submodule(NULL,"QuasiNewton","opt_module");
- * Matrix data;
- * size_t dim_param_space;
- *
- * ...
- * QuasiNewton opt;
- * opt.Init(obj_function, data, dim_param_space, opt_module);
- * ...
- * opt.Eval(init_pt);
- * // init_pt contains the optimal point found
- * @endcode
- *
- */
-class QuasiNewton {
- private:
- size_t dimension_;
- arma::mat data_;
- long double (*func_ptr_)(const arma::vec&, const arma::mat&, arma::vec&);
-
- public:
- QuasiNewton(){ }
- ~QuasiNewton(){ }
-
- void Init(long double (*fun)(const arma::vec&, const arma::mat&, arma::vec&),
- arma::mat& data) {
- data_ = data;
- func_ptr_ = fun;
- dimension_ = mlpack::CLI::GetParam<int>("opt/param_space_dim");
- }
-
- const arma::mat& data() {
- return data_;
- }
-
- size_t dimension() {
- return dimension_;
- }
-
- void Eval(arma::vec& pt);
- void LineSearch_(arma::vec& pold, long double fold, arma::vec& grad,
- arma::vec& xi, arma::vec& pnew, long double& f_min,
- long double maximum_step_length);
-};
-
-#endif
Copied: mlpack/trunk/src/mlpack/methods/mog/optimizers.hpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/optimizers.h)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/optimizers.hpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/optimizers.hpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,153 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file optimizers.hpp
+ *
+ * Declares classes for two types of optimizer
+ *
+ */
+#ifndef __MLPACK_METHODS_MOG_OPTIMIZERS_HPP
+#define __MLPACK_METHODS_MOG_OPTIMIZERS_HPP
+
+#include <mlpack/core.h>
+
+PARAM_STRING("method", "The method used to optimize", "opt", "");
+
+PARAM_INT_REQ("param_space_dim", "The dimension of the parameter space.", "opt");
+PARAM_INT("MAX_FUNC_EVAL", "The maximum number of function evaluations\
+ allowed to the NelderMead optimizer (defaults to 50000)", "opt", 50000);
+
+PARAM_INT("func_evals", "The number of function evaluations taken by the algorithm", "opt", 0);
+PARAM_INT("MAX_ITERS", "The maximum number of iterations allowed to the function", "opt", 200);
+PARAM_INT("iters", "The number of iterations the algorithm actually went through", "opt", 0);
+
+PARAM_DOUBLE("EPSILON", "Value of epsilon.", "opt", 3.0e-8);
+PARAM_DOUBLE("TOLERANCE", "Tolerance for the minimum movement for the parameter value.", "opt", 1.0e-5);
+PARAM_DOUBLE("gtol", "Tolerance value for the gradient of the function", "opt", 1.0e-7);
+PARAM_DOUBLE("MAX_STEP_SIZE", "The maximum step size in the direction of the gradient.", "opt", 100.0);
+PARAM_DOUBLE("tolerance", "Undocumented parameter", "opt", 1.0e-5);
+
+PARAM_MODULE("opt", "This file contains two optimizers.");
+
+namespace mlpack {
+namespace gmm {
+
+/**
+ * An optimizer using the Nelder Mead method,
+ * also known as the polytope or the simplex
+ * method.
+ *
+ * It does multivariate minimization of an
+ * objective function. If it is optimizing in
+ * 'd' dimensions, it would require 'd+1'
+ * starting points.
+ *
+ * Example use:
+ *
+ * @code
+ * double init_pts[d+1][d];
+ * size_t number_of_function_evaluations;
+ * struct datanode *opt_module = fx_submodule(NULL,"NelderMead","opt_module");
+ * Matrix data;
+ * size_t dim_param_space;
+ *
+ * ...
+ * NelderMead opt;
+ * opt.Init(obj_function, data, dim_param_space, opt_module);
+ * ...
+ * opt.Eval(init_pts);
+ * // init_pts[0] contaings the optimal point found
+ * @endcode
+ *
+ */
+class NelderMead {
+ private:
+ size_t dimension_;
+ arma::mat data_;
+ long double (*func_ptr_)(const arma::vec&, const arma::mat&);
+
+ public:
+ NelderMead() { }
+ ~NelderMead() { }
+
+ void Init(long double (*fun)(const arma::vec&, const arma::mat&),
+ arma::mat& data) {
+ data_ = data;
+ func_ptr_ = fun;
+ dimension_ = mlpack::CLI::GetParam<int>("opt/param_space_dim");
+ }
+
+ const arma::mat& data() {
+ return data_;
+ }
+
+ size_t dimension() {
+ return dimension_;
+ }
+
+ void Eval(arma::mat& pts);
+ long double ModSimplex_(arma::mat& pts, arma::vec& y,
+ arma::vec& psum, size_t ihi, float fac);
+};
+
+/**
+ * An optimizer using the Quasi Newton method,
+ * also known as the variable metrics
+ * method.
+ *
+ * It does multivariate minimization of an
+ * objective function using only the function
+ * value and the gradients.
+ *
+ * Example use:
+ *
+ * @code
+ * double init_pt[d];
+ * size_t number_of_iters;
+ * struct datanode *opt_module = fx_submodule(NULL,"QuasiNewton","opt_module");
+ * Matrix data;
+ * size_t dim_param_space;
+ *
+ * ...
+ * QuasiNewton opt;
+ * opt.Init(obj_function, data, dim_param_space, opt_module);
+ * ...
+ * opt.Eval(init_pt);
+ * // init_pt contains the optimal point found
+ * @endcode
+ *
+ */
+class QuasiNewton {
+ private:
+ size_t dimension_;
+ arma::mat data_;
+ long double (*func_ptr_)(const arma::vec&, const arma::mat&, arma::vec&);
+
+ public:
+ QuasiNewton(){ }
+ ~QuasiNewton(){ }
+
+ void Init(long double (*fun)(const arma::vec&, const arma::mat&, arma::vec&),
+ arma::mat& data) {
+ data_ = data;
+ func_ptr_ = fun;
+ dimension_ = mlpack::CLI::GetParam<int>("opt/param_space_dim");
+ }
+
+ const arma::mat& data() {
+ return data_;
+ }
+
+ size_t dimension() {
+ return dimension_;
+ }
+
+ void Eval(arma::vec& pt);
+ void LineSearch_(arma::vec& pold, long double fold, arma::vec& grad,
+ arma::vec& xi, arma::vec& pnew, long double& f_min,
+ long double maximum_step_length);
+};
+
+}; // namespace gmm
+}; // namespace mlpack
+
+#endif // __MLPACK_METHODS_MOG_OPTIMIZERS_HPP
Deleted: mlpack/trunk/src/mlpack/methods/mog/phi.h
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/phi.h 2011-10-31 15:41:27 UTC (rev 10089)
+++ mlpack/trunk/src/mlpack/methods/mog/phi.h 2011-10-31 15:54:31 UTC (rev 10090)
@@ -1,123 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file phi.h
- *
- * This file computes the Gaussian probability
- * density function
- */
-#include <mlpack/core.h>
-
-/**
- * Calculates the multivariate Gaussian probability density function
- *
- * Example use:
- * @code
- * Vector x, mean;
- * Matrix cov;
- * ....
- * long double f = phi(x, mean, cov);
- * @endcode
- */
-long double phi(const arma::vec& x,
- const arma::vec& mean,
- const arma::mat& cov) {
- long double cdet, f;
- double exponent;
- size_t dim;
- arma::mat cinv = inv(cov);
- arma::vec diff, tmp;
-
- dim = x.n_elem;
- cdet = det(cov);
-
- if (cdet < 0)
- cdet = -cdet;
-
- diff = mean - x;
- tmp = cinv * diff;
- exponent = dot(diff, tmp);
-
- long double tmp1, tmp2;
- tmp1 = 1 / pow(2 * M_PI, dim / 2);
- tmp2 = 1 / sqrt(cdet);
- f = (tmp1 * tmp2 * exp(-exponent / 2));
-
- return f;
-}
-
-/**
- * Calculates the univariate Gaussian probability density function
- *
- * Example use:
- * @code
- * double x, mean, var;
- * ....
- * long double f = phi(x, mean, var);
- * @endcode
- */
-long double phi(const double x, const double mean, const double var) {
- return exp(-1.0 * ((x - mean) * (x - mean) / (2 * var)))
- / sqrt(2 * M_PI * var);
-}
-
-/**
- * Calculates the multivariate Gaussian probability density function
- * and also the gradients with respect to the mean and the variance
- *
- * Example use:
- * @code
- * Vector x, mean, g_mean, g_cov;
- * ArrayList<Matrix> d_cov; // the dSigma
- * ....
- * long double f = phi(x, mean, cov, d_cov, &g_mean, &g_cov);
- * @endcode
- */
-long double phi(const arma::vec& x,
- const arma::vec& mean,
- const arma::mat& cov,
- const std::vector<arma::mat>& d_cov,
- arma::vec& g_mean,
- arma::vec& g_cov) {
- long double cdet, f;
- double exponent;
- size_t dim;
- arma::mat cinv = inv(cov);
- arma::vec diff, tmp;
-
- dim = x.n_elem;
- cdet = det(cov);
-
- if (cdet < 0)
- cdet = -cdet;
-
- diff = mean - x;
- tmp = cinv * diff;
- exponent = dot(diff, tmp);
-
- long double tmp1, tmp2;
- tmp1 = 1 / pow(2 * M_PI, dim / 2);
- tmp2 = 1 / sqrt(cdet);
- f = (tmp1 * tmp2 * exp(-exponent / 2));
-
- // Calculating the g_mean values which would be a (1 X dim) vector
- g_mean = f * tmp;
-
- // Calculating the g_cov values which would be a (1 X (dim*(dim+1)/2)) vector
- arma::vec g_cov_tmp(d_cov.size());
- for (size_t i = 0; i < d_cov.size(); i++) {
- arma::vec tmp_d;
- arma::mat inv_d;
- long double tmp_d_cov_d_r;
-
- tmp_d = d_cov[i] * tmp;
- tmp_d_cov_d_r = dot(tmp_d,tmp);
- inv_d = cinv * d_cov[i];
- for (size_t j = 0; j < dim; j++)
- tmp_d_cov_d_r += inv_d(j, j);
- g_cov_tmp[i] = f * tmp_d_cov_d_r / 2;
- }
-
- g_cov = g_cov_tmp;
-
- return f;
-}
Copied: mlpack/trunk/src/mlpack/methods/mog/phi.hpp (from rev 10083, mlpack/trunk/src/mlpack/methods/mog/phi.h)
===================================================================
--- mlpack/trunk/src/mlpack/methods/mog/phi.hpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/mog/phi.hpp 2011-10-31 15:54:31 UTC (rev 10090)
@@ -0,0 +1,129 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file phi.hpp
+ *
+ * This file computes the Gaussian probability
+ * density function
+ */
+#include <mlpack/core.h>
+
+namespace mlpack {
+namespace gmm {
+
+/**
+ * Calculates the multivariate Gaussian probability density function
+ *
+ * Example use:
+ * @code
+ * Vector x, mean;
+ * Matrix cov;
+ * ....
+ * long double f = phi(x, mean, cov);
+ * @endcode
+ */
+long double phi(const arma::vec& x,
+ const arma::vec& mean,
+ const arma::mat& cov) {
+ long double cdet, f;
+ double exponent;
+ size_t dim;
+ arma::mat cinv = inv(cov);
+ arma::vec diff, tmp;
+
+ dim = x.n_elem;
+ cdet = det(cov);
+
+ if (cdet < 0)
+ cdet = -cdet;
+
+ diff = mean - x;
+ tmp = cinv * diff;
+ exponent = dot(diff, tmp);
+
+ long double tmp1, tmp2;
+ tmp1 = 1 / pow(2 * M_PI, dim / 2);
+ tmp2 = 1 / sqrt(cdet);
+ f = (tmp1 * tmp2 * exp(-exponent / 2));
+
+ return f;
+}
+
+/**
+ * Calculates the univariate Gaussian probability density function
+ *
+ * Example use:
+ * @code
+ * double x, mean, var;
+ * ....
+ * long double f = phi(x, mean, var);
+ * @endcode
+ */
+long double phi(const double x, const double mean, const double var) {
+ return exp(-1.0 * ((x - mean) * (x - mean) / (2 * var)))
+ / sqrt(2 * M_PI * var);
+}
+
+/**
+ * Calculates the multivariate Gaussian probability density function
+ * and also the gradients with respect to the mean and the variance
+ *
+ * Example use:
+ * @code
+ * Vector x, mean, g_mean, g_cov;
+ * ArrayList<Matrix> d_cov; // the dSigma
+ * ....
+ * long double f = phi(x, mean, cov, d_cov, &g_mean, &g_cov);
+ * @endcode
+ */
+long double phi(const arma::vec& x,
+ const arma::vec& mean,
+ const arma::mat& cov,
+ const std::vector<arma::mat>& d_cov,
+ arma::vec& g_mean,
+ arma::vec& g_cov) {
+ long double cdet, f;
+ double exponent;
+ size_t dim;
+ arma::mat cinv = inv(cov);
+ arma::vec diff, tmp;
+
+ dim = x.n_elem;
+ cdet = det(cov);
+
+ if (cdet < 0)
+ cdet = -cdet;
+
+ diff = mean - x;
+ tmp = cinv * diff;
+ exponent = dot(diff, tmp);
+
+ long double tmp1, tmp2;
+ tmp1 = 1 / pow(2 * M_PI, dim / 2);
+ tmp2 = 1 / sqrt(cdet);
+ f = (tmp1 * tmp2 * exp(-exponent / 2));
+
+ // Calculating the g_mean values which would be a (1 X dim) vector
+ g_mean = f * tmp;
+
+ // Calculating the g_cov values which would be a (1 X (dim*(dim+1)/2)) vector
+ arma::vec g_cov_tmp(d_cov.size());
+ for (size_t i = 0; i < d_cov.size(); i++) {
+ arma::vec tmp_d;
+ arma::mat inv_d;
+ long double tmp_d_cov_d_r;
+
+ tmp_d = d_cov[i] * tmp;
+ tmp_d_cov_d_r = dot(tmp_d,tmp);
+ inv_d = cinv * d_cov[i];
+ for (size_t j = 0; j < dim; j++)
+ tmp_d_cov_d_r += inv_d(j, j);
+ g_cov_tmp[i] = f * tmp_d_cov_d_r / 2;
+ }
+
+ g_cov = g_cov_tmp;
+
+ return f;
+}
+
+}; // namespace gmm
+}; // namespace mlpack
More information about the mlpack-svn
mailing list