From 126de015791458346cd58aff9bcc20220eb9d805 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 28 Mar 2024 17:16:27 +0100 Subject: Initial commit. --- .gitmodules | 6 +++ least_squares.cpp | 119 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ pcg-cpp | 1 + randutils | 1 + 4 files changed, 127 insertions(+) create mode 100644 .gitmodules create mode 100644 least_squares.cpp create mode 160000 pcg-cpp create mode 160000 randutils diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..a14394b --- /dev/null +++ b/.gitmodules @@ -0,0 +1,6 @@ +[submodule "randutils"] + path = randutils + url = https://gist.github.com/imneme/540829265469e673d045 +[submodule "pcg-cpp"] + path = pcg-cpp + url = https://github.com/imneme/pcg-cpp diff --git a/least_squares.cpp b/least_squares.cpp new file mode 100644 index 0000000..8e7d085 --- /dev/null +++ b/least_squares.cpp @@ -0,0 +1,119 @@ +#include +#include +#include + +#include "pcg-cpp/include/pcg_random.hpp" +#include "randutils/randutils.hpp" + + +template +using Vector = Eigen::Matrix; + +template +using Matrix = Eigen::Matrix; + +template +class Model { +private: + Matrix A; + Vector b; + +public: + template + Model(Real σ, unsigned N, unsigned M, Generator& r) : A(M, N), b(M) { + std::normal_distribution aDistribution(0, 1); + + for (unsigned i = 0; i < M; i++) { + for (unsigned j =0; j < N; j++) { + A(i, j) = aDistribution(r); + } + } + + std::normal_distribution bDistribution(0, σ); + + for (unsigned i = 0; i < M; i++) { + b(i) = bDistribution(r); + } + } + + const unsigned N() { + return A.cols(); + } + + const unsigned M() { + return A.rows(); + } + + const Vector V(const Vector& x) { + return A * x + b; + } + + const Matrix dV(const Vector& x) { + return A; + } + +// const Real ddV(const Vector& x) { +// return Matrix::Zero(; +// } + + const Real H(const Vector& x) { + return V(x).squaredNorm(); + } + + const Vector dH(const Vector& x) { + return dV(x).transpose() * V(x); + } + + const Matrix ddH(const Vector& x) { + return dV(x).transpose() * dV(x); + } + + const Vector ∇H(const Vector& x){ + return dH(x) - dH(x).dot(x) * x / x.squaredNorm(); + } + + const Matrix HessH(const Vector& x) { + Matrix hess = ddH(x) - x.dot(dH(x)) * Matrix::Identity(N(), N()); + return hess - (hess * x) * x.transpose() / x.squaredNorm(); + } +}; + +using Rng = randutils::random_generator; +using Real = double; + +int main(int argc, char* argv[]) { + unsigned N = 10; + unsigned M = 10; + Real σ = 1; + + int opt; + + while ((opt = getopt(argc, argv, "N:M:s:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned)atof(optarg); + break; + case 'M': + M = (unsigned)atof(optarg); + break; + case 's': + σ = atof(optarg); + break; + default: + exit(1); + } + } + + Rng r; + + Model leastSquares(σ, N, M, r.engine()); + + Vector x = Vector::Zero(N); + x(0) = N; + + std::cout << leastSquares.H(x) << std::endl; + std::cout << leastSquares.∇H(x) << std::endl; + std::cout << leastSquares.HessH(x) << std::endl; + + return 0; +} diff --git a/pcg-cpp b/pcg-cpp new file mode 160000 index 0000000..428802d --- /dev/null +++ b/pcg-cpp @@ -0,0 +1 @@ +Subproject commit 428802d1a5634f96bcd0705fab379ff0113bcf13 diff --git a/randutils b/randutils new file mode 160000 index 0000000..5afb443 --- /dev/null +++ b/randutils @@ -0,0 +1 @@ +Subproject commit 5afb4439c23a6c060eda72121bff8bf9da59591a -- cgit v1.2.3-70-g09d2