Use Eigen's SVD to solve linear least squares

_code/CMakeLists.txt

@@ -16,3 +16,4 @@ add_subdirectory(notes)
 add_subdirectory(pset_03)
 add_subdirectory(pset_04)
 add_subdirectory(pset_05/cpp)
+add_subdirectory(pset_07/cpp)

_code/pset_07/cpp/CMakeLists.txt

+#add_library(sim
+#    glut_grapher.cpp
+#    glut_grapher.h
+#)
+#target_link_libraries(sim PUBLIC common ${OPENGL_LIBRARIES} ${GLUT_LIBRARY})
+#target_include_directories(sim PUBLIC
+#    ${OPENGL_INCLUDE_DIRS}
+#    ${GLUT_INCLUDE_DIRS}
+#)
+
+add_executable(svd
+    svd.cpp
+)
+target_link_libraries(svd shared_settings shared_code)
+
+#add_executable(diffusion
+#    diffusion.cpp
+#)
+#target_link_libraries(diffusion sim)

_code/pset_07/cpp/svd.cpp

+#include "xorshift.h"
+#include <Eigen/Dense>
+#include <iostream>
+#include <random>
+#include <vector>
+
+using namespace std;
+using namespace Eigen;
+
+//--------------------------------------------------------------------------------------------------
+XorShift64 rg1;
+std::random_device rd{};
+std::mt19937 rg2{rd()};
+
+//--------------------------------------------------------------------------------------------------
+Eigen::ArrayXd draw_n_uniform(uint32_t n) {
+    Eigen::ArrayXd result(n);
+    for (uint32_t i = 0; i < n; ++i) {
+        result[i] = rg1.draw_double();
+    }
+    return result;
+}
+
+//--------------------------------------------------------------------------------------------------
+// TODO Box-Muller transform my uniform samples instead of using STL
+Eigen::ArrayXd draw_n_normal(uint32_t n, double std_deviation) {
+    std::normal_distribution<> normal{0, std_deviation};
+    Eigen::ArrayXd result(n);
+    for (uint32_t i = 0; i < n; ++i) {
+        result[i] = normal(rg2);
+    }
+    return result;
+}
+
+//--------------------------------------------------------------------------------------------------
+int main() {
+    uint32_t constexpr n_warmup = 1000;
+    for (uint32_t i = 0; i < n_warmup; ++i) {
+        rg1.advance();
+    }
+
+    uint32_t constexpr n_samples = 100;
+    double const std_deviation = 0.5;
+    Eigen::ArrayXd samples = draw_n_uniform(n_samples);
+    Eigen::ArrayXd errors = draw_n_normal(n_samples, std_deviation);
+    Eigen::ArrayXd values = 2 + 3 * samples + errors;
+
+    Eigen::MatrixXd mat(n_samples, 2);
+    for (uint32_t i = 0; i < n_samples; ++i) {
+        mat(i, 0) = samples[i];
+        mat(i, 1) = 1;
+    }
+
+    // Need at least ThinU and ThinV to use mat_svd.solve.
+    JacobiSVD<MatrixXd> mat_svd(mat, ComputeThinU | ComputeThinV);
+    std::cout << "singular values:\n" << mat_svd.singularValues() << '\n';
+    //cout << "left singular vectors:\n" << mat_svd.matrixU() << '\n';
+    //cout << "right singular vectors:\n" << mat_svd.matrixV() << '\n';
+    Eigen::VectorXd solution = mat_svd.solve(values.matrix());
+    std::cout << "solution:\n" << solution << '\n';
+}