Added quadratic, linear, constant approx. for functions, added newtonian 2nd order method for optimization

MLPP/LinReg/LinReg.cpp

@@ -34,6 +34,36 @@ namespace MLPP{
         return Evaluate(x);
     }
 
+    void LinReg::NewtonRaphson(double learning_rate, int max_epoch, bool UI){
+        LinAlg alg;
+        Reg regularization;
+        double cost_prev = 0;
+        int epoch = 1;
+        forwardPass();   
+        while(true){
+            cost_prev = Cost(y_hat, outputSet);
+                
+            std::vector<double> error = alg.subtraction(y_hat, outputSet);
+
+            // Calculating the weight gradients (2nd derivative)
+            std::vector<double> first_derivative = alg.mat_vec_mult(alg.transpose(inputSet), error);
+            std::vector<std::vector<double>> second_derivative = alg.matmult(alg.transpose(inputSet), inputSet);
+            weights = alg.subtraction(weights, alg.scalarMultiply(learning_rate/n, alg.mat_vec_mult(alg.transpose(alg.inverse(second_derivative)), first_derivative)));
+            weights = regularization.regWeights(weights, lambda, alpha, reg);
+ 
+            // Calculating the bias gradients (2nd derivative)
+            bias -= learning_rate * alg.sum_elements(error) / n; // We keep this the same. The 2nd derivative is just [1].
+            forwardPass();
+                
+            if(UI) { 
+                Utilities::CostInfo(epoch, cost_prev, Cost(y_hat, outputSet));
+                Utilities::UI(weights, bias); 
+            }
+            epoch++;
+            if(epoch > max_epoch) { break; }
+        }
+    }
+
     void LinReg::gradientDescent(double learning_rate, int max_epoch, bool UI){
         LinAlg alg;
         Reg regularization;
@@ -59,7 +89,6 @@ namespace MLPP{
                 Utilities::UI(weights, bias); 
             }
             epoch++;
-            
             if(epoch > max_epoch) { break; }
         }
     }

MLPP/LinReg/LinReg.hpp

@@ -17,6 +17,7 @@ namespace MLPP{
             LinReg(std::vector<std::vector<double>> inputSet, std::vector<double> outputSet, std::string reg = "None", double lambda = 0.5, double alpha = 0.5);
             std::vector<double> modelSetTest(std::vector<std::vector<double>> X);
             double modelTest(std::vector<double> x);
+            void NewtonRaphson(double learning_rate, int max_epoch, bool UI);
             void gradientDescent(double learning_rate, int max_epoch, bool UI = 1);
             void SGD(double learning_rate, int max_epoch, bool UI = 1);
             void MBGD(double learning_rate, int max_epoch, int mini_batch_size, bool UI = 1);

MLPP/NumericalAnalysis/NumericalAnalysis.cpp

@@ -5,6 +5,7 @@
 //
 
 #include "NumericalAnalysis.hpp"
+#include "LinAlg/LinAlg.hpp"
 #include <iostream>
 
 namespace MLPP{
@@ -20,6 +21,18 @@ namespace MLPP{
         return (function(x + eps) -  2 * function(x) + function(x - eps)) / (eps * eps);
     }
 
+    double  NumericalAnalysis::constantApproximation(double(*function)(double), double c){
+        return function(c);
+    }
+
+    double  NumericalAnalysis::linearApproximation(double(*function)(double), double c, double x){
+        return constantApproximation(function, c) + numDiff(function, c) * (x - c);
+    }
+
+    double NumericalAnalysis::quadraticApproximation(double(*function)(double), double c, double x){
+        return linearApproximation(function, c, x) + 0.5 * numDiff_2(function, c) * (x - c) * (x - c);
+    }
+
     double NumericalAnalysis::numDiff(double(*function)(std::vector<double>), std::vector<double> x, int axis){
         // For multivariable function analysis. 
         // This will be used for calculating Jacobian vectors. 
@@ -77,4 +90,18 @@ namespace MLPP{
         }
         return hessian;
     }
+
+    double NumericalAnalysis::constantApproximation(double(*function)(std::vector<double>), std::vector<double> c){
+        return function(c);
+    }
+
+    double NumericalAnalysis::linearApproximation(double(*function)(std::vector<double>), std::vector<double> c, std::vector<double> x){
+        LinAlg alg;
+        return constantApproximation(function, c) + alg.matmult(alg.transpose({jacobian(function, c)}), {alg.subtraction(x, c)})[0][0];
+    }
+
+    double NumericalAnalysis::quadraticApproximation(double(*function)(std::vector<double>), std::vector<double> c, std::vector<double> x){
+        LinAlg alg;
+        return linearApproximation(function, c, x) + 0.5 * alg.matmult({(alg.subtraction(x, c))}, alg.matmult(hessian(function, c), alg.transpose({alg.subtraction(x, c)})))[0][0];
+    }
 }

MLPP/NumericalAnalysis/NumericalAnalysis.hpp

@@ -18,6 +18,10 @@ namespace MLPP{
             double numDiff(double(*function)(double), double x);
             double numDiff_2(double(*function)(double), double x);
 
+            double constantApproximation(double(*function)(double), double c);
+            double linearApproximation(double(*function)(double), double c, double x);
+            double quadraticApproximation(double(*function)(double), double c, double x);
+
             double numDiff(double(*function)(std::vector<double>), std::vector<double> x, int axis);
             double numDiff_2(double(*function)(std::vector<double>), std::vector<double> x, int axis1, int axis2);
 
@@ -25,6 +29,11 @@ namespace MLPP{
 
             std::vector<double> jacobian(double(*function)(std::vector<double>), std::vector<double> x); // Indeed, for functions with scalar outputs the Jacobians will be vectors.
             std::vector<std::vector<double>> hessian(double(*function)(std::vector<double>), std::vector<double> x);
+
+            double constantApproximation(double(*function)(std::vector<double>), std::vector<double> c);
+            double linearApproximation(double(*function)(std::vector<double>), std::vector<double> c, std::vector<double> x);
+            double quadraticApproximation(double(*function)(std::vector<double>), std::vector<double> c, std::vector<double> x);
+
         
     };
 }

README.md

@@ -149,7 +149,11 @@ The result will be the model's predictions for the entire dataset.
         - Multivariate Functions 
     2. Jacobian Vector Calculator
     3. Hessian Matrix Calculator
-    3. Newton-Raphson Method
+    4. Function approximator
+        - Constant Approximation
+        - Linear Approximation 
+        - Quadratic Approximation
+    5. Newton-Raphson Method
 14. ***Linear Algebra Module***
 15. ***Statistics Module***
 16. ***Data Processing Module***

main.cpp

@@ -49,8 +49,12 @@
 using namespace MLPP;
 
 
+// double f(double x){
+//     return x*x*x + 2*x - 2; 
+// }
+
 double f(double x){
-    return x*x*x + 2*x - 2; 
+    return cos(x);
 }
 /*
     y = x^3 + 2x - 2
@@ -77,7 +81,7 @@ double f_mv(std::vector<double> x){
     Where x, y = x[0], x[1], this function is defined as:
     f(x, y) = x^4 + xy^3 + yz^2
 
-    ∂f/∂x = 3x^2 + 3y^2 
+    ∂f/∂x = 4x^3 + 3y^2 
     ∂^2f/∂x^2 = 6x 
 
     ∂f/∂z = 2zy
@@ -178,17 +182,37 @@ int main() {
     // UniLinReg model(inputSet, outputSet);
     // alg.printVector(model.modelSetTest(inputSet));
 
-    // MULIVARIATE LINEAR REGRESSION
+    // // MULIVARIATE LINEAR REGRESSION
     // std::vector<std::vector<double>> inputSet = {{1,2,3,4,5,6,7,8,9,10}, {3,5,9,12,15,18,21,24,27,30}};
     // std::vector<double> outputSet = {2,4,6,8,10,12,14,16,18,20};
+
     // LinReg model(alg.transpose(inputSet), outputSet); // Can use Lasso, Ridge, ElasticNet Reg
-    // model.normalEquation(); 
+
     // model.gradientDescent(0.001, 30000, 1);
     // model.SGD(0.001, 30000, 1);
     // model.MBGD(0.001, 10000, 2, 1);
+    // model.normalEquation(); 
+
     // alg.printVector(model.modelSetTest((alg.transpose(inputSet))));
     // std::cout << "ACCURACY: " << 100 * model.score() << "%" << std::endl;
 
+
+    // std::cout << "Total epoch num: 300" << std::endl;
+    // std::cout << "Method: 1st Order w/ Jacobians" << std::endl;
+
+    // LinReg model(alg.transpose(inputSet), outputSet); // Can use Lasso, Ridge, ElasticNet Reg
+
+    // model.gradientDescent(0.001, 300, 0);
+
+
+    // std::cout << "--------------------------------------------" << std::endl;
+    // std::cout << "Total epoch num: 300" << std::endl;
+    // std::cout << "Method: Newtonian 2nd Order w/ Hessians" << std::endl;
+    // LinReg model2(alg.transpose(inputSet), outputSet); 
+
+    // model2.NewtonRaphson(1.5, 300, 0);
+
+
     // // LOGISTIC REGRESSION
     // std::vector<std::vector<double>> inputSet; 
     // std::vector<double> outputSet; 
@@ -474,8 +498,6 @@ int main() {
     // data.getImage("../../Data/apple.jpeg", chicken);
     // alg.printVector(chicken);
 
-    // // TESTING QR DECOMP. EXAMPLE VIA WIKIPEDIA. SEE https://en.wikipedia.org/wiki/QR_decomposition.
-
     // std::vector<std::vector<double>> P = {{12, -51, 4}, {6, 167, -68}, {-4, 24, -41}};
     // alg.printMatrix(P);
 
@@ -504,6 +526,10 @@ int main() {
     // Checks for numerical analysis class.
     NumericalAnalysis numAn;
 
+    //std::cout << numAn.quadraticApproximation(f, 0, 1) << std::endl;
+
+    std::cout << numAn.quadraticApproximation(f_mv, {0, 0, 0}, {1, 1, 1}) << std::endl;
+
     // std::cout << numAn.numDiff(&f, 1) << std::endl;
     // std::cout << numAn.newtonRaphsonMethod(&f, 1, 1000) << std::endl;
 
@@ -513,9 +539,9 @@ int main() {
 
     //std::cout << numAn.numDiff_2(&f, 2) << std::endl;
 
-    std::cout << numAn.numDiff_2(&f_mv, {2, 2, 500}, 2, 2) << std::endl;
-    std::cout << "Our Hessian." << std::endl;
-    alg.printMatrix(numAn.hessian(&f_mv, {2, 2, 500}));
+    // std::cout << numAn.numDiff_2(&f_mv, {2, 2, 500}, 2, 2) << std::endl;
+    // std::cout << "Our Hessian." << std::endl;
+    // alg.printMatrix(numAn.hessian(&f_mv, {2, 2, 500}));
 
     return 0;
 }