diff options
| -rw-r--r-- | fits.cpp | 119 |
1 files changed, 61 insertions, 58 deletions
@@ -1,6 +1,6 @@ #include <getopt.h> #include <iomanip> -#include<list> +#include <list> #include "pcg-cpp/include/pcg_random.hpp" #include "randutils/randutils.hpp" @@ -12,80 +12,73 @@ using Rng = randutils::random_generator<pcg32>; using Real = double; using Vector = Eigen::Matrix<Real, Eigen::Dynamic, 1>; -class fitModel { -public: - unsigned N; +/* +inline Real basisFunctions(unsigned i, Real x) const { + if (i == 0) { + return 1; + } else { + return pow(x, i); + } +} +*/ - fitModel(unsigned N) : N(N) {} +inline Real basisFunctions(unsigned N, unsigned i, Real x) { + return abs(x - i / (N - 1.0)); +} - /* - virtual Real basisFunctions(unsigned i, Real x) const { - if (i == 0) { - return 1; - } else { - return pow(x, i); - } - } - */ +Real value(const Vector& coeffs, Real x) { + Real v = 0; - Real basisFunctions(unsigned i, Real x) const { - return abs(x - i / (N - 1.0)); + for (unsigned j = 0; j < coeffs.size(); j++) { + v += coeffs[j] * basisFunctions(coeffs.size(), j, x); } - Real value(const Vector& coeffs, Real x) const { - Real v = 0; + return v; +} - for (unsigned j = 0; j < N; j++) { - v += coeffs[j] * basisFunctions(j, x); - } +Real cost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) { + Real c = 0; - return v; + for (const std::tuple<Real, Real>& xy : data) { + Real x = std::get<0>(xy); + Real y = std::get<1>(xy); + c += 0.5 * pow(value(coeffs, x) - y, 2); } - Real cost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) const { - Real c = 0; + return c; +} - for (std::tuple<Real, Real> xy : data) { +Vector dCost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) { + Vector dC = Vector::Zero(coeffs.size()); + + for (unsigned j = 0; j < coeffs.size(); j++) { + for (const std::tuple<Real, Real>& xy : data) { Real x = std::get<0>(xy); Real y = std::get<1>(xy); - c += 0.5 * pow(value(coeffs, x) - y, 2); + dC[j] += (value(coeffs, x) - y) * basisFunctions(coeffs.size(), j, x); } - - return c; } - Vector dCost(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs) const { - Vector dC = Vector::Zero(N); - - for (unsigned j = 0; j < N; j++) { - for (std::tuple<Real, Real> xy : data) { - Real x = std::get<0>(xy); - Real y = std::get<1>(xy); - dC[j] += (value(coeffs, x) - y) * basisFunctions(j, x); - } - } - - return dC; - } -}; + return dC; +} -Vector gradientDescent(const fitModel& M, const std::list<std::tuple<Real, Real>>& data, const Vector& a₀, unsigned maxSteps, Real ε = 1e-12) { +Vector gradientDescent(const std::list<std::tuple<Real, Real>>& data, const Vector& a₀, unsigned maxSteps, Real ε = 1e-12) { Vector xₜ = a₀; - Real Hₜ = M.cost(data, a₀); + Real Hₜ = cost(data, a₀); Real α = 1.0; Real m; Vector ∇H; unsigned steps = 0; while ( - ∇H = M.dCost(data, xₜ), m = ∇H.squaredNorm(), + ∇H = dCost(data, xₜ), m = ∇H.squaredNorm(), m > ε && steps < maxSteps ) { Vector xₜ₊₁; Real Hₜ₊₁; while ( - xₜ₊₁ = xₜ - α * ∇H, Hₜ₊₁ = M.cost(data, xₜ₊₁), + xₜ₊₁ = xₜ - α * ∇H, Hₜ₊₁ = cost(data, xₜ₊₁), Hₜ₊₁ > Hₜ - 0.5 * α * m ) { α /= 2; @@ -100,7 +93,23 @@ Vector gradientDescent(const fitModel& M, const std::list<std::tuple<Real, Real> return xₜ; } -Vector stochasticGradientDescent(const fitModel& M, const std::list<std::tuple<Real, Real>>& data, const Vector& a₀, Rng& r, unsigned maxSteps, Real ε = 1e-12) { +Vector dCostRand(const std::list<std::tuple<Real, Real>>& data, const Vector& coeffs, Real batchP, Rng& r) { + Vector dC = Vector::Zero(coeffs.size()); + + for (unsigned j = 0; j < coeffs.size(); j++) { + for (const std::tuple<Real, Real>& xy : data) { + if (batchP > r.uniform(0.0, 1.0)) { + Real x = std::get<0>(xy); + Real y = std::get<1>(xy); + dC[j] += (value(coeffs, x) - y) * basisFunctions(coeffs.size(), j, x); + } + } + } + + return dC; +} + +Vector stochasticGradientDescent(std::list<std::tuple<Real, Real>>& data, const Vector& a₀, Rng& r, unsigned maxSteps, Real ε = 1e-12) { Vector xₜ = a₀; Real α = 1e-3; Real m; @@ -108,14 +117,10 @@ Vector stochasticGradientDescent(const fitModel& M, const std::list<std::tuple<R unsigned steps = 0; while ( - ∇H = M.dCost(data, xₜ), m = ∇H.squaredNorm(), + ∇H = dCostRand(data, xₜ, 0.1, r), m = ∇H.squaredNorm(), m > ε && steps < maxSteps ) { - Vector γ(M.N); - for (Real& γi : γ) { - γi = r.variate<Real, std::normal_distribution>(0, 1e-4); - } - xₜ = xₜ - α * ∇H + γ; + xₜ = xₜ - α * ∇H; steps++; } @@ -137,7 +142,7 @@ int main(int argc, char* argv[]) { unsigned N = 10; unsigned M = 10; Real σ = 0.2; - unsigned maxSteps = 1e5; + unsigned maxSteps = 1e8; int opt; @@ -162,7 +167,7 @@ int main(int argc, char* argv[]) { Rng r; - std::list<std::tuple<Real, Real>> data = generateData([](Real x) {return sin(2 * M_PI * x);}, M, σ, r); + std::list<std::tuple<Real, Real>> data = generateData([](Real x) {return x;}, M, σ, r); std::cout << std::setprecision(15); @@ -171,10 +176,8 @@ int main(int argc, char* argv[]) { } std::cout << std::endl; - fitModel model(N); - Vector a₀ = Vector::Zero(N); - Vector a = gradientDescent(model, data, a₀, maxSteps); + Vector a = stochasticGradientDescent(data, a₀, r, maxSteps); for (unsigned i = 0; i < N; i++) { std::cout << a[i] << " "; |