diff --git a/bin/abc_lab b/bin/abc_lab
index 1a58c21..45555df 100755
Binary files a/bin/abc_lab and b/bin/abc_lab differ
diff --git a/include/modules/neural_networks/activation_functions/Activation_ReLU.h b/include/modules/neural_networks/activation_functions/Activation_ReLU.h
index add7906..2427503 100644
--- a/include/modules/neural_networks/activation_functions/Activation_ReLU.h
+++ b/include/modules/neural_networks/activation_functions/Activation_ReLU.h
@@ -11,12 +11,29 @@ namespace neural_networks{
 	template <typename T>
 	struct Activation_ReLU{
 
-		//utils::Matrix<T> inputs;
+		utils::Matrix<T> _inputs;
 		utils::Matrix<T> outputs;
 
+		utils::Matrix<T> dinputs;
+
         void forward(const utils::Matrix<T>& inputs){
+        	_inputs = inputs;
 		    outputs = numerics::matclip_low(inputs, T{0});
         }
+        void backward(const utils::Matrix<T>& dvalues){
+        	// Since we need to modify the original variable,
+        	// let's make a copy of the values first
+        	dinputs = dvalues;
+
+        	// Zero gradients where input values were negative
+        	for (uint64_t i = 0; i < dinputs.rows(); ++i){
+        		for (uint64_t j = 0; j < dinputs.cols(); ++j){
+        			if (_inputs(i,j) <= T{0}){
+        				dinputs(i,j) = T{0};
+        			}
+        		}
+        	}
+        }
 	};
 
 
diff --git a/include/modules/neural_networks/layers/Dense_Layer.h b/include/modules/neural_networks/layers/Dense_Layer.h
index 4acc229..ecf9faa 100644
--- a/include/modules/neural_networks/layers/Dense_Layer.h
+++ b/include/modules/neural_networks/layers/Dense_Layer.h
@@ -12,11 +12,15 @@ namespace neural_networks{
 	template <typename T>
 	struct Dense_Layer{
 
-			//utils::Matrix<T> _inputs;
+			utils::Matrix<T> _inputs;
 			utils::Matrix<T> weights;
 			utils::Vector<T> biases;
 			utils::Matrix<T> outputs;
 
+			utils::Matrix<T> dweights;
+			utils::Vector<T> dbiases;
+			utils::Matrix<T> dinputs;
+
 			
 			// Default Constructor
 			Dense_Layer() = default;
@@ -29,9 +33,20 @@ namespace neural_networks{
 	        	
 	        }
 
-	        void forward(utils::Matrix<T>& inputs){
+	        void forward(const utils::Matrix<T>& inputs){
+	        	_inputs = inputs;
 			    outputs = numerics::matadd(numerics::matmul_auto(inputs, weights), biases, "row");
 	        }
+
+	        void backward(const utils::Matrix<T>& dvalues){
+	        	// Gradients on parameters
+	        	dweights  = numerics::matmul(numerics::transpose(_inputs), dvalues);
+	        	dbiases = numerics::matsum(dvalues, "row");
+	        	//Gradient on values
+	        	dinputs = numerics::matmul(dvalues, numerics::transpose(weights));
+
+	        }
+
 	};
 
 
diff --git a/include/numerics/matrandom.h b/include/numerics/matrandom.h
new file mode 100644
index 0000000..b64093d
--- /dev/null
+++ b/include/numerics/matrandom.h
@@ -0,0 +1,55 @@
+#ifndef _matrandom_n_
+#define _matrandom_n_
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+
+namespace numerics{
+	
+	template <typename T>
+	void inplace_matrandom_add(utils::Matrix<T>& A, const T lower, const T higher) {
+	    
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+	    utils::Matrix<T> B;
+	    B.random(rows,cols, lower, higher);
+
+	    numerics::inplace_matadd_mat(A, B);
+	}
+
+	template <typename T>
+	void inplace_matrandom_mul(utils::Matrix<T>& A, const T lower, const T higher) {
+	    
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+	    utils::Matrix<T> B;
+	    B.random(rows,cols, lower, higher);
+
+	    for (uint64_t i = 0; i < rows; ++i){
+	    	for (uint64_t j = 0; j < cols; ++j){
+	    		A(i,j) *= B(i,j);
+	    	}
+	    }
+	}
+
+
+	template <typename T>
+	utils::Matrix<T> matrandom_add(const utils::Matrix<T>& A, const T lower, const T higher) {
+	    
+	    const uint64_t rows = A.rows();
+	    const uint64_t cols = A.cols();
+
+	    utils::Matrix<T> B = A;
+
+	    numerics::inplace_matadd_mat(B, lower, higher);
+
+	    return B;
+	}
+
+
+} // namespace numerics
+
+#endif // _matrandom_n_
\ No newline at end of file
diff --git a/include/numerics/numerics.h b/include/numerics/numerics.h
index 27f2a10..c3135e4 100644
--- a/include/numerics/numerics.h
+++ b/include/numerics/numerics.h
@@ -8,6 +8,7 @@
 #include "./numerics/vecmax.h"
 #include "./numerics/veclog.h"
 #include "./numerics/vecargmax.h"
+#include "./numerics/vecrandom.h"
 #include "./numerics/initializers/eye.h"
 #include "./numerics/matequal.h"
 #include "./numerics/transpose.h"
@@ -18,6 +19,7 @@
 #include "./numerics/matdiv.h"
 #include "./numerics/matvec.h"
 #include "./numerics/matadd.h"
+#include "./numerics/matrandom.h"
 #include "./numerics/matsubtract.h"
 #include "./numerics/matsum.h"
 #include "./numerics/matclip.h"
diff --git a/include/numerics/vecrandom.h b/include/numerics/vecrandom.h
new file mode 100644
index 0000000..75e60b5
--- /dev/null
+++ b/include/numerics/vecrandom.h
@@ -0,0 +1,48 @@
+#ifndef _vecrandom_n_
+#define _vecrandom_n_
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+#include "./core/omp_config.h"
+
+namespace numerics{
+	
+	template <typename T>
+	void inplace_vecrandom_add(utils::Vector<T>& a, const T lower, const T higher) {
+	    
+	    const uint64_t N = a.size();
+
+	    utils::Vector<T> b;
+	    b.random(N, lower, higher);
+
+	    a.inplace_add(b);
+	}
+
+	template <typename T>
+	void inplace_vecrandom_mul(utils::Vector<T>& a, const T lower, const T higher) {
+	    
+	    const uint64_t N = a.size();
+
+	    utils::Vector<T> b;
+	    b.random(N, lower, higher);
+
+	    a.inplace_multiply(b);
+
+	}
+
+
+	template <typename T>
+	utils::Vector<T> vecrandom_add(const utils::Vector<T>& a, const T lower, const T higher) {
+	    
+	    const uint64_t N = a.size();
+
+	    utils::Vector<T> b;
+	    inplace_vecrandom_add(b, lower, higher);
+
+	    return b;
+	}
+
+
+} // namespace numerics
+
+#endif // _vecrandom_n_
\ No newline at end of file
diff --git a/include/utils/matrix.h b/include/utils/matrix.h
index 8929e72..7a41f76 100644
--- a/include/utils/matrix.h
+++ b/include/utils/matrix.h
@@ -23,7 +23,6 @@ namespace utils{
 	class Matrix{
 public:
 	    Matrix() : rows_(0), cols_(0), data_() {} // Default constructor
-#include <random>
 	   // Constructor to initialize matrix with rows × cols and a fill value
 	    Matrix(uint64_t rows, uint64_t cols, const T& value = T()) 
 	    : rows_(rows), cols_(cols), data_(rows * cols, value) {}
diff --git a/obj/main.d b/obj/main.d
index b7ed3f9..7fffdb3 100644
--- a/obj/main.d
+++ b/obj/main.d
@@ -7,19 +7,20 @@ obj/main.o: src/main.cpp include/./core/omp_config.h \
  include/./numerics/exp.h include/./numerics/log.h \
  include/./numerics/vecclip.h include/./numerics/vecexp.h \
  include/./numerics/vecmax.h include/./numerics/veclog.h \
- include/./numerics/vecargmax.h include/./numerics/initializers/eye.h \
- include/./numerics/matequal.h include/./numerics/abs.h \
- include/./numerics/transpose.h include/./numerics/inverse.h \
+ include/./numerics/vecargmax.h include/./numerics/vecrandom.h \
+ include/./numerics/initializers/eye.h include/./numerics/matequal.h \
+ include/./numerics/abs.h include/./numerics/transpose.h \
+ include/./numerics/inverse.h \
  include/./numerics/inverse/inverse_gauss_jordan.h \
  include/./numerics/inverse/inverse_lu.h include/./decomp/lu.h \
  include/./numerics/matmul.h include/./numerics/matscalar.h \
  include/./numerics/matmax.h include/./numerics/matdiv.h \
  include/./numerics/matvec.h include/./numerics/matadd.h \
- include/./numerics/matsubtract.h include/./numerics/matsum.h \
- include/./numerics/matclip.h include/./numerics/matexp.h \
- include/./numerics/matlog.h include/./numerics/matdot.h \
- include/./numerics/matargmax.h include/./numerics/min.h \
- include/./numerics/interpolation1d.h \
+ include/./numerics/matrandom.h include/./numerics/matsubtract.h \
+ include/./numerics/matsum.h include/./numerics/matclip.h \
+ include/./numerics/matexp.h include/./numerics/matlog.h \
+ include/./numerics/matdot.h include/./numerics/matargmax.h \
+ include/./numerics/min.h include/./numerics/interpolation1d.h \
  include/./numerics/interpolation1d/interpolation1d_barycentric.h \
  include/./numerics/interpolation1d/interpolation1d_base.h \
  include/./numerics/interpolation1d/interpolation1d_cubic_spline.h \
@@ -57,6 +58,7 @@ include/./numerics/vecexp.h:
 include/./numerics/vecmax.h:
 include/./numerics/veclog.h:
 include/./numerics/vecargmax.h:
+include/./numerics/vecrandom.h:
 include/./numerics/initializers/eye.h:
 include/./numerics/matequal.h:
 include/./numerics/abs.h:
@@ -71,6 +73,7 @@ include/./numerics/matmax.h:
 include/./numerics/matdiv.h:
 include/./numerics/matvec.h:
 include/./numerics/matadd.h:
+include/./numerics/matrandom.h:
 include/./numerics/matsubtract.h:
 include/./numerics/matsum.h:
 include/./numerics/matclip.h:
diff --git a/obj/main.o b/obj/main.o
index daee6c6..8511dde 100644
Binary files a/obj/main.o and b/obj/main.o differ
diff --git a/src/main.cpp b/src/main.cpp
index 1ac4389..e85ed5f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -44,7 +44,7 @@ int main(int argc, char const *argv[])
     utils::Vf best_dense_1_biases = dense1.biases;
     utils::Mf best_dense_2_weights = dense2.weights;
     utils::Vf best_dense_2_biases = dense2.biases;
-    utils::Mf matrRND;
+
     utils::Vf vectRND;
 
     utils::Vector<uint64_t> predections;
@@ -53,20 +53,15 @@ int main(int argc, char const *argv[])
 
 
 
-        for (uint64_t i = 0; i < 10000; ++i){
+        for (uint64_t i = 0; i < 1000; ++i){
 
 
             // Generate a new set of weights for iteration
-            matrRND.random(dense1.weights.rows(), dense1.weights.cols(), -1, 1);
-            numerics::inplace_matscalar(matrRND,0.05f);
-            numerics::inplace_matadd(dense1.weights, matrRND);
+            numerics::inplace_matrandom_mul(dense1.weights,0.98f, 1.02f);
+            numerics::inplace_vecrandom_mul(dense1.biases,0.98f, 1.02f);
 
-            dense1.biases.random(dense1.biases.size(), -1, 1);
-
-            matrRND.random(dense2.weights.rows(), dense2.weights.cols(), -1, 1);
-            numerics::inplace_matscalar(matrRND,0.05f);
-            dense2.weights.random(dense2.weights.rows(), dense2.weights.cols(), -1, 1);
-            dense2.biases.random(dense2.biases.size(), -1, 1);
+            numerics::inplace_matrandom_mul(dense2.weights,0.98f, 1.02f);
+            numerics::inplace_vecrandom_mul(dense2.biases,0.98f, 1.02f);
 
             // Perform a forward pass of the training data through this layer
             dense1.forward(X);
@@ -95,36 +90,29 @@ int main(int argc, char const *argv[])
                 best_dense_2_weights = dense2.weights;
                 best_dense_2_biases = dense2.biases;
                 lowest_loss = loss;
+            } else{
+                //std::cout << "HERE" << std::endl;
+                dense1.weights = best_dense_1_weights;
+                dense1.biases = best_dense_1_biases;
+                dense2.weights = best_dense_2_weights;
+                dense2.biases = best_dense_2_biases;
             }
 
 
         }
 
 
-
-
-
-
-
-
-    for (int i = 0; i < 5; ++i){
-        std::cout << activation2.outputs.get_row(i) << std::endl;
-    }
-
     std::cout << loss << std::endl;
-
-    
-    //predections.print();
-    
-    
-
-
-
-    //predections.print();
-    //class_targets.print();
-
     std::cout << accuracy << std::endl;
 
+    utils::Mf dvalues{{1, 1, 1},
+                    {2, 2, 2},
+                    {3, 3, 3}};
+    utils::Vf biases{2,3,0.5};
+
+    utils::Vf dbiases = numerics::matsum(dvalues, "row");
+
+    dbiases.print();