From 5e0d141ffef12a4623a1e33cc11ce01d900c1ce5 Mon Sep 17 00:00:00 2001
From: Marcin Kirsz <mkirsz@ed.ac.uk>
Date: Tue, 18 Feb 2025 21:22:51 +0000
Subject: [PATCH] Memory manager

---
 include/tadah/models/Algorithm.h              |  30 +
 include/tadah/models/ea.h                     | 122 +++-
 include/tadah/models/ekm.h                    |   9 +-
 include/tadah/models/linear_regressor.h       | 129 ++--
 .../models/memory/IModelsWorkspaceManager.h   |  66 ++
 .../models/memory/ModelsWorkspaceManager.h    |  78 +++
 include/tadah/models/memory/OLSWorkspace.h    |  89 +++
 include/tadah/models/memory/SVDWorkspace.h    |  70 ++
 include/tadah/models/ols.h                    | 514 +++-----------
 include/tadah/models/ols_impl.h               | 311 +++++++++
 include/tadah/models/ridge_regression.h       | 110 ---
 include/tadah/models/svd.h                    | 322 +++------
 include/tadah/models/svd_impl.h               |  79 +++
 src/ModelsWorkspaceManager.cpp                |  54 ++
 src/OLSWorkspace.cpp                          | 167 +++++
 src/SVDWorkspace.cpp                          |  76 +++
 src/ols.cpp                                   | 165 -----
 tests/test_ols.cpp                            | 636 +++++++++---------
 tests/test_ridge_regression.cpp               |  31 -
 tests/test_svd.cpp                            | 270 ++++++--
 20 files changed, 1883 insertions(+), 1445 deletions(-)
 create mode 100644 include/tadah/models/Algorithm.h
 create mode 100644 include/tadah/models/memory/IModelsWorkspaceManager.h
 create mode 100644 include/tadah/models/memory/ModelsWorkspaceManager.h
 create mode 100644 include/tadah/models/memory/OLSWorkspace.h
 create mode 100644 include/tadah/models/memory/SVDWorkspace.h
 create mode 100644 include/tadah/models/ols_impl.h
 delete mode 100644 include/tadah/models/ridge_regression.h
 create mode 100644 include/tadah/models/svd_impl.h
 create mode 100644 src/ModelsWorkspaceManager.cpp
 create mode 100644 src/OLSWorkspace.cpp
 create mode 100644 src/SVDWorkspace.cpp
 delete mode 100644 tests/test_ridge_regression.cpp

diff --git a/include/tadah/models/Algorithm.h b/include/tadah/models/Algorithm.h
new file mode 100644
index 0000000..d03c5f7
--- /dev/null
+++ b/include/tadah/models/Algorithm.h
@@ -0,0 +1,30 @@
+#ifndef TADAH_MODELS_ALGORITHM_H
+#define TADAH_MODELS_ALGORITHM_H
+
+namespace tadah {
+namespace models {
+
+/**
+ * @enum Algorithm
+ * @brief Enum to select the algorithm used for solving computational problems.
+ *
+ * This enum is used by various classes (e.g., OLS) to specify the algorithm to be used.
+ *
+ * **Algorithms:**
+ * - **GELS (`DGELS`):** Uses QR or LQ factorization, fast but less robust for ill-conditioned systems; suitable for overdetermined systems.
+ * - **GELSD (`DGELSD`):** Uses Singular Value Decomposition (SVD) approach, handling rank-deficient and ill-conditioned systems; suitable for both overdetermined and underdetermined systems.
+ * - **Pseudoinverse:** Custom implementation using SVD to compute the pseudoinverse; allows for explicit control over singular values and regularization; suitable for underdetermined systems.
+ * - **NormalEquations:** Solves the normal equations \f$ A^\top A x = A^\top b \f$ using Cholesky decomposition.
+ */
+enum class Algorithm {
+    GELS,           ///< Uses DGELS
+    GELSD,          ///< Uses DGELSD
+    Pseudoinverse,  ///< Uses custom SVD-based pseudoinverse implementation
+    NormalEquations ///< Uses Cholesky decomposition to solve normal equations
+};
+
+} // namespace models
+} // namespace tadah
+
+#endif // TADAH_MODELS_ALGORITHM_H
+
diff --git a/include/tadah/models/ea.h b/include/tadah/models/ea.h
index 2f6987d..fd19b3f 100644
--- a/include/tadah/models/ea.h
+++ b/include/tadah/models/ea.h
@@ -1,9 +1,86 @@
+#ifndef SVD_MULTIPLIER_H
+#define SVD_MULTIPLIER_H
+#include <tadah/models/svd.h>
+#include <tadah/core/lapack.h>
+#include <vector>
+#include <algorithm>
+
+/**
+ * @class SVDMultiplier
+ * @brief Class to compute A * x using the SVD decomposition of A without allocating memory during computation.
+ */
+class SVDMultiplier {
+public:
+    /**
+     * @brief Constructor that initializes the multiplier with the SVD decomposition.
+     *
+     * @param svd The SVD object containing the decomposition of A.
+     */
+    explicit SVDMultiplier(const tadah::models::SVD& svd)
+        : m_(svd.rows()), n_(svd.cols()), min_mn_(std::min(m_, n_)),
+          U_(svd.getU()), S_(svd.getS()), VT_(svd.getVT()),
+          y_(min_mn_), z_(min_mn_) {}
+
+    /**
+     * @brief Computes the matrix-vector product A * x using the SVD decomposition.
+     *
+     * @param x Pointer to the input vector x (of size n).
+     * @param result Pointer to the output vector where A * x will be stored (of size m).
+     */
+    void compute(const double* x, double* result) {
+        // Step 1: Compute y = V^T * x
+        {
+            char trans = 'N'; // V^T is already transposed
+            int m = min_mn_;
+            int n = n_;
+            double alpha = 1.0;
+            int lda = m;
+            int incx = 1;
+            double beta = 0.0;
+            int incy = 1;
+
+            dgemv_(&trans, &m, &n, &alpha, VT_, &lda, x, &incx, &beta, y_.data(), &incy);
+        }
+
+        // Step 2: Compute z = S * y (element-wise multiplication)
+        for (int i = 0; i < min_mn_; ++i) {
+            z_[i] = S_[i] * y_[i];
+        }
+
+        // Step 3: Compute result = U * z
+        {
+            char trans = 'N';
+            int m = m_;
+            int n = min_mn_;
+            double alpha = 1.0;
+            int lda = m;
+            int incx = 1;
+            double beta = 0.0;
+            int incy = 1;
+
+            dgemv_(&trans, &m, &n, &alpha, U_, &lda, z_.data(), &incx, &beta, result, &incy);
+        }
+    }
+
+private:
+    int m_;                  // Number of rows in A
+    int n_;                  // Number of columns in A
+    int min_mn_;             // Minimum of m and n
+
+    const double* U_;        // Pointer to U matrix
+    const double* S_;        // Pointer to singular values
+    const double* VT_;       // Pointer to V^T matrix
+
+    std::vector<double> y_;  // Intermediate vector y (size min_mn_)
+    std::vector<double> z_;  // Intermediate vector z (size min_mn_)
+};
+#endif
+
 #ifndef M_EA_H
 #define M_EA_H
 
 #include <tadah/core/config.h>
 #include <tadah/models/svd.h>
-#include <tadah/models/ridge_regression.h>
 #include <tadah/models/ols.h>
 #include <tadah/core/core_types.h>
 #include <tadah/core/lapack.h>
@@ -11,11 +88,11 @@
 class EA {
   private:
     const Config &config;
-    const SVD &svd;
+    tadah::models::SVD &svd;
     t_type &T;
     int verbose;
   public:
-    EA(const Config &c, const SVD &svd, t_type &T):
+    EA(const Config &c, tadah::models::SVD &svd, t_type &T):
       config(c),
       svd(svd),
       T(T),
@@ -24,7 +101,8 @@ class EA {
 
     int run(double &alpha, double &beta) {
       // Phi = USV^T
-      size_t N=svd.shapeA().first;
+      size_t N=std::max(svd.rows(), svd.cols());
+      size_t sizeS = std::min(svd.rows(), svd.cols());
       size_t maxsteps=40;
       const double EPS2 = 5e-12;
       int convergence_status = 1;
@@ -39,39 +117,46 @@ class EA {
 
       // square of the S matrix contains all eigenvalues
       // of Phi^T*Phi = (VSU^T)(USV^T) = VS^2V^T = S^2
-      double *eval = new double[svd.sizeS()];
-      for (size_t i=0; i<svd.sizeS(); ++i)
+      double *eval = new double[sizeS];
+      for (size_t i=0; i<sizeS; ++i)
         eval[i] = s[i]*s[i];
 
       if (verbose) printf("%9s  %9s  %9s  %9s  %9s  %9s\n","del alpha", "del beta", "alpha", "beta", "lambda", "gamma");
       // D is a filter factors matrix
-      double *d = new double[svd.sizeS()];
+      double *d = new double[sizeS];
 
       int outprec=config.get<int>("OUTPREC");
 
-      double *t = new double[svd.shapeA().first];
-      double *x = new double[svd.shapeA().second];
-      t_type m_n((size_t)svd.sizeS());
-      OLS::Workspace ws;
-      auto shapeA = svd.shapeA();
-      int m = shapeA.first;
-      int n = shapeA.second;
+      t_type m_n(sizeS);
+      int m = svd.rows();
+      int n = svd.cols();
+      t_type Tpred(m);
     double rcond = config.size("OALGO")==2 ? config.get<double>("OALGO",1) : 1e-8;
+    //
+        // Create a workspace manager
+        tadah::models::memory::ModelsWorkspaceManager workspaceManager;
+
+        // Create OLS instance with the workspace manager
+        tadah::models::OLS ols(workspaceManager);
+
+        SVDMultiplier svdMultiplier(svd);
       while (test1>EPS2 || test2>EPS2) {
 
         // regularised least square problem (Tikhonov Regularization):
         /*RidgeRegression::solve(svd,T,m_n,lambda,rcond);*/
-        OLS::solve_with_svd_with_workspace(m, n, T, m_n, svd, ws, rcond, lambda);
+        /*OLS::solve_with_svd_with_workspace(m, n, T, m_n, svd, ws, rcond, lambda);*/
+        ols.solveWithPseudoinverse(m, n, T, m_n, svd.getSVDWorkspace(), rcond, lambda);
 
         double gamma = 0.0;
-        for (size_t j=0; j<svd.sizeS(); j++) {
+        for (size_t j=0; j<sizeS; j++) {
           if (beta*eval[j] > std::numeric_limits<double>::min()) {
             gamma += (beta*eval[j])/(beta*eval[j]+alpha);
           }
         }
         alpha = gamma / (m_n*m_n);
 
-        t_type Tpred = svd.predict(m_n,t,x);
+        /*t_type Tpred = svd.predict(m_n,t,x);*/
+        svdMultiplier.compute(m_n.data(), Tpred.data());
         double temp = (T-Tpred)*(T-Tpred);
 
         beta = (static_cast<double>(N)-gamma)/temp;
@@ -95,9 +180,6 @@ class EA {
       delete [] eval;
       delete [] d;
 
-      delete[] t;
-      delete[] x;
-
       return convergence_status;
     };
 };
diff --git a/include/tadah/models/ekm.h b/include/tadah/models/ekm.h
index fa4540f..c150ad0 100644
--- a/include/tadah/models/ekm.h
+++ b/include/tadah/models/ekm.h
@@ -61,12 +61,13 @@ class EKM {
             }
             eps *= max_val;
 
-            SVD svd(K);
+            tadah::models::SVD svd;
+            svd.compute(K);
 
             double *s = svd.getS();
             double *vt = svd.getVT();
-            t_type S(s,svd.sizeS());
-            Matrix VT(vt,svd.shapeVT().first, svd.shapeVT().second);
+            t_type S(s,svd.cols());
+            Matrix VT(vt,svd.cols(), svd.cols());
 
             // First count number of non zero singular values, then
             // filter out corresponding right singular vectors from
@@ -74,7 +75,7 @@ class EKM {
             size_t nonzero_sv=0;
             for (size_t i=0; i<S.size(); ++i) if (s[i]>eps) nonzero_sv++;
 
-            EKM_mat = Matrix(nonzero_sv, svd.shapeVT().second);
+            EKM_mat = Matrix(nonzero_sv, svd.cols());
 
             for (size_t i=0; i<S.size(); ++i) {
                 if (s[i]>eps) {
diff --git a/include/tadah/models/linear_regressor.h b/include/tadah/models/linear_regressor.h
index 54df74a..9ea2c14 100644
--- a/include/tadah/models/linear_regressor.h
+++ b/include/tadah/models/linear_regressor.h
@@ -6,25 +6,26 @@
 #include <tadah/core/config.h>
 #include <tadah/core/core_types.h>
 #include <tadah/models/ea.h>
-#include <tadah/models/svd.h>
 #include <tadah/models/ols.h>
-#include <tadah/models/ridge_regression.h>
+#include <tadah/models/svd.h>
 
 #include <vector>
 
 /**
  * @class LinearRegressor
  * @brief Performs linear regression using Singular Value Decomposition (SVD).
- * 
- * Utilizes SVD for training, supporting both ordinary least squares and ridge regression.
+ *
+ * Utilizes SVD for training, supporting both ordinary least squares and ridge
+ * regression.
  */
 class LinearRegressor {
 
   /**
    * @brief Train the model using the provided configuration and data matrices.
-   * 
-   * This method supports both ordinary least squares and ridge regression based on the lambda value.
-   * 
+   *
+   * This method supports both ordinary least squares and ridge regression based
+   * on the lambda value.
+   *
    * @param config Configuration object containing training parameters.
    * @param Phi Design matrix.
    * @param T Target matrix.
@@ -32,35 +33,37 @@ class LinearRegressor {
    * @param Sigma Matrix to store the covariance matrix.
    */
 public:
-  static void train(Config &config, phi_type &Phi, t_type &T,
-                    t_type &weights, Matrix &Sigma) {
+  static void train(Config &config, phi_type &Phi, t_type &T, t_type &weights,
+                    Matrix &Sigma) {
 
     int verbose = config.get<int>("VERBOSE");
-    double lambda = config.get<double>("LAMBDA",0);
-    int oalgo =  config.get<int>("OALGO",0);
-    double rcond = config.size("OALGO")==2 ? config.get<double>("OALGO",1) : 1e-8;
+    double lambda = config.get<double>("LAMBDA", 0);
+    int oalgo = config.get<int>("OALGO", 0);
+    double rcond =
+        config.size("OALGO") == 2 ? config.get<double>("OALGO", 1) : 1e-8;
     weights.resize(Phi.cols());
 
-    if (std::abs(lambda+1) < std::numeric_limits<double>::min() && oalgo != 3) {
+    if (std::abs(lambda + 1) < std::numeric_limits<double>::min() &&
+        oalgo != 3) {
       throw std::runtime_error("LAMBDA -1 requires OALGO 3");
     }
 
-    if (lambda > 0 && (oalgo==1 || oalgo==2)) {
+    if (lambda > 0 && oalgo != 3 ) {
       // Augment the design matrix with the regularization term
-      size_t j=0;
-      for (size_t i=Phi.rows()-Phi.cols(); i<Phi.rows();++i) {
-        Phi(i,j++) = std::sqrt(lambda);
+      size_t j = 0;
+      for (size_t i = Phi.rows() - Phi.cols(); i < Phi.rows(); ++i) {
+        Phi(i, j++) = std::sqrt(lambda);
       }
     }
 
-    if (oalgo==1) {
-      OLS::solve(Phi, T, weights, rcond, OLS::Algorithm::GELSD);
-    }
-    else if (oalgo==2) {
-      OLS::solve(Phi, T, weights, rcond, OLS::Algorithm::GELS);
-    }
-    else if (oalgo==3) {
-      SVD svd = SVD(Phi);;
+    tadah::models::OLS ols;
+    if (oalgo == 1) {
+      ols.solve(Phi, T, weights, rcond, tadah::models::Algorithm::GELSD);
+    } else if (oalgo == 2) {
+      ols.solve(Phi, T, weights, rcond, tadah::models::Algorithm::GELS);
+    } else if (oalgo == 3) {
+      tadah::models::SVD svd;
+      svd.compute(Phi);
       if (lambda < 0) {
         double alpha = config.get<double>("ALPHA");
         double beta = config.get<double>("BETA");
@@ -76,30 +79,32 @@ public:
       Sigma = get_sigma(svd, lambda);
       config_add_sigma(config, Sigma);
 
-      if (verbose) std::cout << std::endl << "REG LAMBDA: " << lambda << std::endl;
-      OLS::solve_with_svd(Phi.rows(), Phi.cols(), T, weights, svd, rcond, lambda);
-    }
-    else if (oalgo==4) {
-      OLS::solve(Phi, T, weights, rcond, OLS::Algorithm::Cholesky);
-    }
-    else {
+      if (verbose)
+        std::cout << std::endl << "REG LAMBDA: " << lambda << std::endl;
+
+      ols.solveWithPseudoinverse(Phi.rows(), Phi.cols(), T, weights,
+                                 svd.getSVDWorkspace(), rcond, lambda);
+    } else if (oalgo == 4) {
+      ols.solve(Phi, T, weights, rcond,
+                tadah::models::Algorithm::NormalEquations);
+    } else {
       throw std::runtime_error("Unsupported OALGO: " + std::to_string(oalgo));
     }
   }
 
-  /** 
-     * @brief Compute and return the covariance matrix.
-     * 
-     * Uses the SVD to calculate the covariance matrix with regularization.
-     * 
-     * @param svd SVD object of the design matrix.
-     * @param lambda Regularization parameter.
-     * @return Covariance matrix.
-     */
-  static Matrix get_sigma(SVD &svd, double lambda) {
+  /**
+   * @brief Compute and return the covariance matrix.
+   *
+   * Uses the SVD to calculate the covariance matrix with regularization.
+   *
+   * @param svd SVD object of the design matrix.
+   * @param lambda Regularization parameter.
+   * @return Covariance matrix.
+   */
+  static Matrix get_sigma(tadah::models::SVD &svd, double lambda) {
     double *VT = svd.getVT();
     double *S = svd.getS();
-    int n = static_cast<int>(svd.shapeVT().first);
+    int n = static_cast<int>(svd.cols());
 
     Matrix Sigma((size_t)n, (size_t)n);
     Sigma.set_zero();
@@ -109,7 +114,8 @@ public:
       double ridge = 1.0 / (S[i] * S[i] + lambda);
       for (int j = 0; j < n; ++j) {
         for (int k = 0; k < n; ++k) {
-          sigma[j + k * n] += VT[j + i * n] * ridge * VT[k + i * n]; // Column-major
+          sigma[j + k * n] +=
+              VT[j + i * n] * ridge * VT[k + i * n]; // Column-major
         }
       }
     }
@@ -117,13 +123,13 @@ public:
   }
 
   /**
-     * @brief Add the covariance matrix to the configuration.
-     * 
-     * Converts the Sigma matrix to a vector and stores it in the configuration.
-     * 
-     * @param config Configuration object.
-     * @param Sigma Covariance matrix to be added.
-     */
+   * @brief Add the covariance matrix to the configuration.
+   *
+   * Converts the Sigma matrix to a vector and stores it in the configuration.
+   *
+   * @param config Configuration object.
+   * @param Sigma Covariance matrix to be added.
+   */
   static void config_add_sigma(Config &config, Matrix &Sigma) {
     t_type Sigma_as_vector(Sigma.data(), Sigma.cols() * Sigma.rows());
     config.add("SIGMA", Sigma.rows());
@@ -131,21 +137,22 @@ public:
   }
 
   /**
-     * @brief Read and reconstruct the covariance matrix from the configuration.
-     * 
-     * Retrieves the Sigma matrix stored in the configuration.
-     * 
-     * @param config Configuration object.
-     * @param Sigma Matrix to store the reconstructed Sigma.
-     * @throws std::runtime_error if Sigma is not computed.
-     */
+   * @brief Read and reconstruct the covariance matrix from the configuration.
+   *
+   * Retrieves the Sigma matrix stored in the configuration.
+   *
+   * @param config Configuration object.
+   * @param Sigma Matrix to store the reconstructed Sigma.
+   * @throws std::runtime_error if Sigma is not computed.
+   */
   static void read_sigma(Config &config, Matrix &Sigma) {
     using V = std::vector<double>;
     size_t N;
     try {
       N = config.get<size_t>("SIGMA");
-    } catch (const std::runtime_error& error) {
-      throw std::runtime_error("\nSigma matrix is not computed.\nHint: It is only computed for LAMBDA != 0.\n");
+    } catch (const std::runtime_error &error) {
+      throw std::runtime_error("\nSigma matrix is not computed.\nHint: It is "
+                               "only computed for LAMBDA != 0.\n");
     }
 
     V Sigma_as_vector(N * N + 1);
diff --git a/include/tadah/models/memory/IModelsWorkspaceManager.h b/include/tadah/models/memory/IModelsWorkspaceManager.h
new file mode 100644
index 0000000..763791f
--- /dev/null
+++ b/include/tadah/models/memory/IModelsWorkspaceManager.h
@@ -0,0 +1,66 @@
+#ifndef TADAH_MODELS_MEMORY_IMODELSWORKSPACEMANAGER_H
+#define TADAH_MODELS_MEMORY_IMODELSWORKSPACEMANAGER_H
+
+#include <tadah/core/memory/IWorkspaceManager.h>
+#include <tadah/models/Algorithm.h> // Include the Algorithm enum
+
+namespace tadah {
+namespace models {
+namespace memory {
+
+class OLSWorkspace; ///< Forward declaration of OLSWorkspace
+class SVDWorkspace; ///< Forward declaration of SVDWorkspace
+
+/**
+ * @interface IModelsWorkspaceManager
+ * @brief Interface for managing workspaces specific to the MODELS module.
+ *
+ * Extends the core IWorkspaceManager interface and provides methods for
+ * obtaining and releasing OLS and SVD workspaces.
+ */
+class IModelsWorkspaceManager : public tadah::core::memory::IWorkspaceManager {
+public:
+    /**
+     * @brief Virtual destructor.
+     */
+    virtual ~IModelsWorkspaceManager() = default;
+
+    /**
+     * @brief Get an OLS workspace suitable for the given problem size and algorithm.
+     *
+     * @param m Number of rows in the matrix.
+     * @param n Number of columns in the matrix.
+     * @param algorithm The algorithm to be used.
+     * @return Pointer to an allocated OLSWorkspace.
+     */
+    virtual OLSWorkspace* getOLSWorkspace(int m, int n, tadah::models::Algorithm algorithm) = 0;
+
+    /**
+     * @brief Release an OLS workspace.
+     *
+     * @param workspace Pointer to the OLSWorkspace to be released.
+     */
+    virtual void releaseOLSWorkspace(OLSWorkspace* workspace) = 0;
+
+    /**
+     * @brief Get an SVD workspace suitable for the given problem size.
+     *
+     * @param m Number of rows in the matrix.
+     * @param n Number of columns in the matrix.
+     * @return Pointer to an allocated SVDWorkspace.
+     */
+    virtual SVDWorkspace* getSVDWorkspace(int m, int n) = 0;
+
+    /**
+     * @brief Release an SVD workspace.
+     *
+     * @param workspace Pointer to the SVDWorkspace to be released.
+     */
+    virtual void releaseSVDWorkspace(SVDWorkspace* workspace) = 0;
+};
+
+} // namespace memory
+} // namespace models
+} // namespace tadah
+
+#endif // TADAH_MODELS_MEMORY_IMODELSWORKSPACEMANAGER_H
diff --git a/include/tadah/models/memory/ModelsWorkspaceManager.h b/include/tadah/models/memory/ModelsWorkspaceManager.h
new file mode 100644
index 0000000..9e38ea4
--- /dev/null
+++ b/include/tadah/models/memory/ModelsWorkspaceManager.h
@@ -0,0 +1,78 @@
+#ifndef TADAH_MODELS_MEMORY_MODELSWORKSPACEMANAGER_H
+#define TADAH_MODELS_MEMORY_MODELSWORKSPACEMANAGER_H
+
+#include <tadah/models/memory/IModelsWorkspaceManager.h>
+#include <tadah/models/Algorithm.h>
+
+namespace tadah {
+namespace models {
+namespace memory {
+
+/**
+ * @class ModelsWorkspaceManager
+ * @brief Implements the IModelsWorkspaceManager interface to manage workspaces.
+ *
+ * Manages OLS and SVD workspaces, providing methods to obtain and release them.
+ */
+class ModelsWorkspaceManager : public IModelsWorkspaceManager {
+public:
+    /**
+     * @brief Constructor.
+     */
+    ModelsWorkspaceManager();
+
+    /**
+     * @brief Destructor.
+     *
+     * Releases all allocated workspaces.
+     */
+    ~ModelsWorkspaceManager() override;
+
+    // OLS workspace methods
+
+    /**
+     * @brief Get an OLS workspace.
+     *
+     * @param m Number of rows in the matrix.
+     * @param n Number of columns in the matrix.
+     * @param algorithm The OLS algorithm to be used.
+     * @return Pointer to an allocated OLSWorkspace.
+     */
+    OLSWorkspace* getOLSWorkspace(int m, int n, tadah::models::Algorithm algorithm) override;
+
+    /**
+     * @brief Release an OLS workspace.
+     *
+     * @param workspace Pointer to the OLSWorkspace to be released.
+     */
+    void releaseOLSWorkspace(OLSWorkspace* workspace) override;
+
+    // SVD workspace methods
+
+    /**
+     * @brief Get an SVD workspace.
+     *
+     * @param m Number of rows in the matrix.
+     * @param n Number of columns in the matrix.
+     * @return Pointer to an allocated SVDWorkspace.
+     */
+    SVDWorkspace* getSVDWorkspace(int m, int n) override;
+
+    /**
+     * @brief Release an SVD workspace.
+     *
+     * @param workspace Pointer to the SVDWorkspace to be released.
+     */
+    void releaseSVDWorkspace(SVDWorkspace* workspace) override;
+
+private:
+    OLSWorkspace* olsWorkspace_; ///< Pointer to the OLSWorkspace.
+    SVDWorkspace* svdWorkspace_; ///< Pointer to the SVDWorkspace.
+};
+
+} // namespace memory
+} // namespace models
+} // namespace tadah
+
+#endif // TADAH_MODELS_MEMORY_MODELSWORKSPACEMANAGER_H
+
diff --git a/include/tadah/models/memory/OLSWorkspace.h b/include/tadah/models/memory/OLSWorkspace.h
new file mode 100644
index 0000000..58bccc0
--- /dev/null
+++ b/include/tadah/models/memory/OLSWorkspace.h
@@ -0,0 +1,89 @@
+#ifndef TADAH_MODELS_MEMORY_OLSWORKSPACE_H
+#define TADAH_MODELS_MEMORY_OLSWORKSPACE_H
+
+#include <tadah/models/Algorithm.h> // For Algorithm
+
+namespace tadah {
+namespace models {
+namespace memory {
+
+/**
+ * @class OLSWorkspace
+ * @brief Manages memory allocations for OLS computations.
+ *
+ * The OLSWorkspace class handles the memory required for various OLS
+ * algorithms. It ensures that memory is allocated and deallocated appropriately
+ * and provides storage for intermediate computations.
+ */
+class OLSWorkspace {
+public:
+  /**
+   * @brief Constructor.
+   */
+  OLSWorkspace();
+
+  /**
+   * @brief Destructor.
+   *
+   * Releases all allocated memory.
+   */
+  ~OLSWorkspace();
+
+  /**
+   * @brief Allocate workspace for the specified problem size and algorithm.
+   *
+   * @param m Number of rows in the matrix.
+   * @param n Number of columns in the matrix.
+   * @param algorithm The OLS algorithm to be used.
+   */
+  void allocate(int m_, int n_, tadah::models::Algorithm algorithm);
+
+  /**
+   * @brief Check if the workspace is sufficient for the given problem size and
+   * algorithm.
+   *
+   * @param m Number of rows in the matrix.
+   * @param n Number of columns in the matrix.
+   * @param algorithm The OLS algorithm to be used.
+   * @return True if the workspace is sufficient, false otherwise.
+   */
+  bool isSufficient(int m_, int n_,
+                    tadah::models::Algorithm algorithm) const;
+
+  // Member variables
+  double *work; ///< Workspace array for computations.
+  int lwork;    ///< Size of the work array.
+  int m;        ///< Number of rows in the matrix.
+  int n;        ///< Number of columns in the matrix.
+  tadah::models::Algorithm
+      algo; ///< Algorithm for which workspace was allocated.
+
+  // Additional variables for specific algorithms
+
+  // For GELSD (DGELSD)
+  double *s;  ///< Singular values array.
+  int *iwork; ///< Integer workspace array.
+
+  // For internal SVD computation in OLS (when Algorithm::SVD is used)
+  double lambda; ///< Regularization parameter.
+  double *U;     ///< U matrix from SVD.
+  double *S;     ///< Singular values from SVD.
+  double *VT;    ///< V^T matrix from SVD.
+  int *iworkSVD; ///< Integer workspace array for SVD.
+  double *UtB;   ///< Intermediate vector for U^T * B.
+
+  // For Cholesky
+  double *AtA; ///< Matrix to store A^T * A.
+  double *Atb; ///< Vector to store A^T * b.
+
+private:
+  // Disable copying
+  OLSWorkspace(const OLSWorkspace &) = delete;
+  OLSWorkspace &operator=(const OLSWorkspace &) = delete;
+};
+
+} // namespace memory
+} // namespace models
+} // namespace tadah
+
+#endif // TADAH_MODELS_MEMORY_OLSWORKSPACE_H
diff --git a/include/tadah/models/memory/SVDWorkspace.h b/include/tadah/models/memory/SVDWorkspace.h
new file mode 100644
index 0000000..8695439
--- /dev/null
+++ b/include/tadah/models/memory/SVDWorkspace.h
@@ -0,0 +1,70 @@
+#ifndef TADAH_MODELS_MEMORY_SVDWORKSPACE_H
+#define TADAH_MODELS_MEMORY_SVDWORKSPACE_H
+
+namespace tadah {
+namespace models {
+namespace memory {
+
+/**
+ * @class SVDWorkspace
+ * @brief Manages memory allocations for SVD computations.
+ *
+ * The SVDWorkspace class handles the memory required for Singular Value Decomposition
+ * computations, ensuring efficient memory usage and reusability.
+ */
+class SVDWorkspace {
+public:
+    /**
+     * @brief Constructor.
+     */
+    SVDWorkspace();
+
+    /**
+     * @brief Destructor.
+     *
+     * Releases all allocated memory.
+     */
+    ~SVDWorkspace();
+
+    /**
+     * @brief Allocate workspace for the specified problem size.
+     *
+     * @param m Number of rows in the matrix.
+     * @param n Number of columns in the matrix.
+     */
+    void allocate(int m_, int n_);
+
+    /**
+     * @brief Check if the workspace is sufficient for the given problem size.
+     *
+     * @param m Number of rows in the matrix.
+     * @param n Number of columns in the matrix.
+     * @return True if the workspace is sufficient, false otherwise.
+     */
+    bool isSufficient(int m_, int n_) const;
+
+    // Member variables for SVD computations
+    int m; ///< Number of rows in the matrix.
+    int n; ///< Number of columns in the matrix.
+
+    double* U;   ///< U matrix from SVD (m x min(m,n)).
+    double* S;   ///< Singular values from SVD (min(m,n)).
+    double* VT;  ///< V^T matrix from SVD (min(m,n) x n).
+
+    // Workspace variables
+    double* work; ///< Workspace array for computations.
+    int lwork;    ///< Size of the work array.
+    int* iwork;   ///< Integer workspace array.
+
+private:
+    // Disable copying
+    SVDWorkspace(const SVDWorkspace&) = delete;
+    SVDWorkspace& operator=(const SVDWorkspace&) = delete;
+};
+
+} // namespace memory
+} // namespace models
+} // namespace tadah
+
+#endif // TADAH_MODELS_MEMORY_SVDWORKSPACE_H
+
diff --git a/include/tadah/models/ols.h b/include/tadah/models/ols.h
index 1c427a9..1971689 100644
--- a/include/tadah/models/ols.h
+++ b/include/tadah/models/ols.h
@@ -1,456 +1,132 @@
-#ifndef OLS_H
-#define OLS_H
+#ifndef TADAH_MODELS_OLS_H
+#define TADAH_MODELS_OLS_H
 
-#include <tadah/core/lapack.h> // Ensure this header includes the necessary LAPACK function declarations
-#include <tadah/core/core_types.h>
+#include <tadah/models/Algorithm.h> // Include the Algorithm enum
+#include <tadah/models/memory/IModelsWorkspaceManager.h>
 #include <cmath>
 #include <stdexcept>
-#include <cstdlib>    // For malloc and free
-#include <algorithm>  // For std::max
-#include <iostream>   // For std::ostream
+#include <cstdlib>
+#include <algorithm>
+#include <iostream>
+
+namespace tadah {
+namespace models {
 
 /**
  * @class OLS
  * @brief Provides functionality for solving Ordinary Least Squares (OLS) problems.
  *
- * Utilizes various algorithms to solve linear systems where the goal is to minimize
- * the Euclidean norm of residuals. Supports the following algorithms:
+ * The OLS class implements various algorithms to solve linear systems where
+ * the goal is to minimize the Euclidean norm of residuals.
  *
  * **Algorithms:**
- * - **GELSD (`DGELSD`):** Uses a singular value decomposition (SVD) approach, handling rank-deficient and ill-conditioned systems.
  * - **GELS (`DGELS`):** Uses QR or LQ factorization, fast but less robust for ill-conditioned systems.
- * - **SVD:** Custom implementation using SVD, allows for explicit control over singular values and regularization.
- * - **Cholesky:** Solves the normal equations \( A^\top A x = A^\top b \) using Cholesky decomposition. Efficient for well-conditioned systems; however, as it squares the condition number, it may be less stable for ill-conditioned systems.
- *
- * **Workspace Management:**
- * - The class provides a `Workspace` struct to manage memory allocations for the computations.
- * - Workspaces are allocated based on the dimensions of the problem and the algorithm selected.
- *
- * **Example:**
- * ```cpp
- * Matrix A; // Initialize with your data
- * Vector B; // Initialize with your data
- * Vector weights(A.cols());
- * OLS::Workspace ws;
- *
- * // Solve using default algorithm (GELSD)
- * OLS::solve(A, B, weights);
- *
- * // Solve using Cholesky algorithm
- * OLS::solve_with_workspace(A, B, weights, ws, -1.0, OLS::Algorithm::Cholesky);
- *
- * // Solve using precomputed SVD and regularization
- * SVD svd(A);
- * double lambda = 0.1;
- * OLS::solve_with_svd(A, B, weights, svd, lambda);
- * ```
+ * - **GELSD (`DGELSD`):** Uses Singular Value Decomposition (SVD) approach, handling rank-deficient and ill-conditioned systems; suitable for both overdetermined and underdetermined systems.
+ * - **Pseudoinverse:** Custom implementation using SVD to compute the pseudoinverse; allows for explicit control over singular values and regularization; suitable for underdetermined systems.
+ * - **NormalEquations:** Solves the normal equations \f$ A^\top A x = A^\top b \f$ using Cholesky decomposition.
  */
 class OLS {
 public:
     /**
-     * @brief Enum to select the algorithm used for solving the OLS problem.
+     * @brief Overload the << operator for the Algorithm enum class.
+     *
+     * @param os The output stream.
+     * @param algo The algorithm to output.
+     * @return Reference to the output stream.
      */
-    enum class Algorithm {
-        GELS,     // Uses DGELS
-        GELSD,    // Uses DGELSD
-        SVD,      // Uses custom SVD implementation
-        Cholesky  // Uses Cholesky decomposition
-    };
-
-    // Overload the << operator for the Algorithm enum class
-    friend std::ostream& operator<<(std::ostream& os, Algorithm algo) {
-        switch (algo) {
-            case Algorithm::GELS:
-                os << "GELS";
-                break;
-            case Algorithm::GELSD:
-                os << "GELSD";
-                break;
-            case Algorithm::SVD:
-                os << "SVD";
-                break;
-            case Algorithm::Cholesky:
-                os << "Cholesky";
-                break;
-            default:
-                os << "Unknown Algorithm";
-                break;
-        }
-        return os;
-    }
+    friend std::ostream& operator<<(std::ostream& os, Algorithm algo);
 
     /**
-     * @brief Workspace struct to hold pre-allocated memory.
+     * @brief Default constructor.
      *
-     * This struct contains the necessary buffers for the computations in different algorithms.
-     * Users can create an instance of `Workspace` and pass it to the `solve_with_workspace` method.
+     * Creates an instance of OLS with its own workspace manager.
      */
-    struct Workspace {
-        // Common workspace variables
-        double* work; // Double workspace array
-        int lwork;    // Size of the work array
-        int m;        // Number of rows in A (for which workspace was allocated)
-        int n;        // Number of columns in A (for which workspace was allocated)
-        Algorithm algo; // Algorithm for which workspace was allocated
-
-        // For DGELSD (GELSD algorithm)
-        double* s;    // Singular values array
-        int* iwork;   // Integer workspace array
-
-        // For SVD algorithm
-        double lambda;
-        double* U;    // U matrix from SVD (m x min(m,n))
-        double* S;    // Singular values array (size min(m,n))
-        double* VT;   // V^T matrix from SVD (min(m,n) x n)
-        int* iwork_svd; // Integer workspace array for SVD
-        double* U_t_b;  // Store result of U^T * b
-
-        // Ownership flags for U, S, VT
-        bool owns_U;
-        bool owns_S;
-        bool owns_VT;
+    OLS();
 
-        // For Cholesky algorithm
-        double* AtA;  // A^T * A matrix (n x n)
-        double* Atb;  // A^T * b vector (n)
-
-        // Constructors and methods (declarations only)
-        Workspace();
-        ~Workspace();
-
-        void allocate(int m_, int n_, Algorithm algorithm = Algorithm::GELSD);
-        void release();
-        bool is_sufficient(int m_, int n_, Algorithm algorithm = Algorithm::GELSD) const;
-    };
-
-    // Public methods
+    /**
+     * @brief Constructor with external workspace manager.
+     *
+     * @param workspaceManager Reference to an external IModelsWorkspaceManager.
+     */
+    explicit OLS(memory::IModelsWorkspaceManager& workspaceManager);
 
-    // Original solve methods
-    template <typename M, typename V>
-    static void solve(M& A, V& B, V& weights);
+    /**
+     * @brief Destructor.
+     *
+     * Deletes the workspace manager if owned.
+     */
+    ~OLS();
 
+    /**
+     * @brief Solve the OLS problem using the default algorithm (GELSD).
+     *
+     * @tparam M Type of the matrix A.
+     * @tparam V Type of the vectors B and weights.
+     * @param A The matrix A.
+     * @param B The vector B.
+     * @param weights Vector to store the solution weights.
+     */
     template <typename M, typename V>
-    static void solve(M& A, V& B, V& weights, double rcond, Algorithm algorithm);
+    void solve(M& A, V& B, V& weights);
 
+    /**
+     * @brief Solve the OLS problem using the specified algorithm.
+     *
+     * @tparam M Type of the matrix A.
+     * @tparam V Type of the vectors B and weights.
+     * @param A The matrix A.
+     * @param B The vector B.
+     * @param weights Vector to store the solution weights.
+     * @param rcond Reciprocal condition number; used to determine effective rank.
+     * @param algorithm The algorithm to use for solving.
+     */
     template <typename M, typename V>
-    static void solve_with_workspace(M& A, V& B, V& weights, Workspace& ws);
+    void solve(M& A, V& B, V& weights, double rcond, Algorithm algorithm);
 
-    template <typename M, typename V>
-    static void solve_with_workspace(M& A, V& B, V& weights, Workspace& ws, double rcond, Algorithm algorithm);
+    /**
+     * @brief Solve the OLS problem using a precomputed SVD (Pseudoinverse method).
+     *
+     * @tparam V1 Type of the vector B.
+     * @tparam V2 Type of the vector weights.
+     * @param m Number of rows in the matrix A.
+     * @param n Number of columns in the matrix A.
+     * @param B The vector B.
+     * @param weights Vector to store the solution weights.
+     * @param svdWorkspace Reference to an SVDWorkspace containing precomputed SVD.
+     * @param rcond Reciprocal condition number; used to determine effective rank (default is -1.0).
+     */
+template <typename V1, typename V2>
+void solveWithPseudoinverse(int m, int n, V1& B, V2& weights, memory::SVDWorkspace& svdWorkspace, double rcond, double lambda);
+template <typename V1, typename V2>
+void solveWithPseudoinverse(int m, int n, V1& B, V2& weights, memory::SVDWorkspace& svdWorkspace);
 
-    // New methods: solve_with_svd
 
-    // Solve using precomputed SVD (without workspace)
-    template <typename V1 ,typename V2, typename SVDType>
-      static void solve_with_svd(int m_, int n_, V1& B, V2& weights, SVDType& svd);
+private:
+    memory::IModelsWorkspaceManager* workspaceManager_; ///< Pointer to the workspace manager.
+    bool ownsWorkspaceManager_; ///< Indicates if this instance owns the workspace manager.
 
-    template <typename V1 ,typename V2, typename SVDType>
-      static void solve_with_svd(int m_, int n_, V1& B, V2& weights, SVDType& svd, double rcond, double lambda);
+    // Internal implementation methods
 
-    // Solve using precomputed SVD with workspace
-    template <typename V1 ,typename V2, typename SVDType>
-      static void solve_with_svd_with_workspace(int m_, int n_, V1& B, V2& weights, SVDType& svd, Workspace& ws);
+    /**
+     * @brief Internal implementation of the solve method.
+     *
+     * @tparam M Type of the matrix A.
+     * @tparam V Type of the vectors B and weights.
+     * @param A The matrix A.
+     * @param B The vector B.
+     * @param weights Vector to store the solution weights.
+     * @param rcond Reciprocal condition number; used to determine effective rank.
+     * @param algorithm The algorithm to use for solving.
+     */
+    template <typename M, typename V>
+    void solveImpl(M& A, V& B, V& weights, double rcond, Algorithm algorithm);
 
-    template <typename V1 ,typename V2, typename SVDType>
-      static void solve_with_svd_with_workspace(int m_, int n_, V1& B, V2& weights, SVDType& svd, Workspace& ws, double rcond, double lambda);
 };
 
-// Implement template methods here (definitions must be in the header)
-
-// Original solve methods
-
-template <typename M, typename V>
-void OLS::solve(M& A, V& B, V& weights) {
-    double rcond = -1.0;
-    Algorithm algorithm = Algorithm::GELSD;
-    solve(A, B, weights, rcond, algorithm);
-}
-
-template <typename M, typename V>
-void OLS::solve(M& A, V& B, V& weights, double rcond, Algorithm algorithm) {
-    Workspace ws;
-    solve_with_workspace(A, B, weights, ws, rcond, algorithm);
-}
-
-template <typename M, typename V>
-void OLS::solve_with_workspace(M& A, V& B, V& weights, Workspace& ws) {
-    double rcond = -1.0;
-    Algorithm algorithm = Algorithm::GELSD;
-    solve_with_workspace(A, B, weights, ws, rcond, algorithm);
-}
-
-// The main 'solve_with_workspace' template method
-template <typename M, typename V>
-void OLS::solve_with_workspace(M& A, V& B, V& weights, Workspace& ws, double rcond, Algorithm algorithm) {
-    int m = A.rows();
-    int n = A.cols();
-
-    // Check for valid algorithm
-    if (algorithm != Algorithm::GELSD && algorithm != Algorithm::GELS &&
-        algorithm != Algorithm::SVD && algorithm != Algorithm::Cholesky) {
-        throw std::invalid_argument("Invalid algorithm specified.");
-    }
-
-    // Check that dimensions of weights match expected size
-    if (weights.size() != static_cast<std::size_t>(n)) {
-        throw std::invalid_argument("weights size is incorrect.");
-    }
-
-    // Check that B has sufficient size
-    int maxmn = std::max(m, n);
-    if (B.size() < static_cast<std::size_t>(maxmn)) {
-        throw std::invalid_argument("B size is insufficient.");
-    }
-
-    // Check if workspace is sufficient; if not, reallocate
-    if (!ws.is_sufficient(m, n, algorithm)) {
-        ws.allocate(m, n, algorithm);
-    }
-
-    int info;
-    if (algorithm == Algorithm::GELSD) {
-        int nrhs = 1;
-        double* a = A.ptr();
-        int lda = m;
-        double* b = B.ptr();
-        int ldb = maxmn;
-        int rank;
-
-        // Solve using DGELSD
-        ::dgelsd_(&m, &n, &nrhs, a, &lda, b, &ldb, ws.s, &rcond, &rank,
-                  ws.work, &ws.lwork, ws.iwork, &info);
-
-        if (info != 0) {
-            throw std::runtime_error("Error in DGELSD: info = " + std::to_string(info));
-        }
-
-        // Copy solution to weights
-        for (int i = 0; i < n; ++i) {
-            weights[i] = b[i];
-        }
-    } else if (algorithm == Algorithm::GELS) {
-        int nrhs = 1;
-        double* a = A.ptr();
-        int lda = m;
-        double* b = B.ptr();
-        int ldb = m;
-        char trans = 'N';
-
-        // Solve using DGELS
-        ::dgels_(&trans, &m, &n, &nrhs, a, &lda, b, &ldb,
-                 ws.work, &ws.lwork, &info);
-
-        if (info != 0) {
-            throw std::runtime_error("Error in DGELS: info = " + std::to_string(info));
-        }
-
-        // Copy solution to weights
-        for (int i = 0; i < n; ++i) {
-            weights[i] = b[i];
-        }
-    } else if (algorithm == Algorithm::SVD) {
-
-        // Custom SVD-based solution
-        int min_mn = std::min(m, n);
-
-        //if (!ws.U || !ws.S || !ws.VT) {
-        if (ws.owns_U || ws.owns_S || ws.owns_VT) {
-            // Perform SVD
-
-            // Ensure workspace is allocated
-            if (!ws.is_sufficient(m, n, algorithm)) {
-                ws.allocate(m, n, algorithm);
-            }
-
-            char jobz = 'S'; // Compute the first min(m,n) singular values and vectors
-            int lda = m;
-            int ldu = m;
-            int ldvt = min_mn; // VT is min(m,n) x n
-
-            // Now perform SVD
-            ::dgesdd_(&jobz, &m, &n, A.ptr(), &lda, ws.S, ws.U, &ldu, ws.VT, &ldvt, ws.work, &ws.lwork, ws.iwork_svd, &info);
-
-            if (info != 0) {
-                throw std::runtime_error("Error in DGESDD: info = " + std::to_string(info));
-            }
-
-        }
-
-        if (ws.U_t_b == nullptr) {
-            ws.U_t_b = new double[min_mn];
-        }
-
-        // Compute Uáµ— * b
-        char trans = 'T'; // Transpose U
-        double alpha = 1.0;
-        double beta = 0.0;
-        int incx = 1;
-        int incy = 1;
-        int lda_u = m;
-        ::dgemv_(&trans, &m, &min_mn, &alpha, ws.U, &lda_u, B.ptr(), &incx, &beta, ws.U_t_b, &incy);
-
-        // Find the maximum singular value
-        double smax = 0.0;
-        for (int i = 0; i < min_mn; ++i) {
-            if (ws.S[i] > smax) {
-                smax = ws.S[i];
-            }
-        }
-
-        // Threshold for singular values
-        double thresh = rcond > 0.0 ? rcond * smax : 0.0;
-
-        // Apply regularization using smax and lambda from workspace
-        for (int i = 0; i < min_mn; ++i) {
-            if (ws.S[i] > thresh) {
-                ws.U_t_b[i] = ws.U_t_b[i] * ws.S[i] / (ws.S[i] * ws.S[i] + ws.lambda);
-            } else {
-                ws.U_t_b[i] = 0.0;
-            }
-        }
-
-        // Compute x = V * (regularized solution)
-        trans = 'T'; // Transpose VT to get V
-        int lda_vt = min_mn;
-        ::dgemv_(&trans, &min_mn, &n, &alpha, ws.VT, &lda_vt, ws.U_t_b, &incx, &beta, weights.ptr(), &incy);
-
-    } else if (algorithm == Algorithm::Cholesky) {
-        // Cholesky-based solution
-
-        // Compute AtA = Aáµ— * A
-        char uplo = 'U'; // Use upper triangle
-        char trans = 'T'; // Transpose A
-        double alpha = 1.0;
-        double beta = 0.0;
-        int lda = m;
-        int ldc = n;
-        ::dsyrk_(&uplo, &trans, &n, &m, &alpha, A.ptr(), &lda, &beta, ws.AtA, &ldc);
-
-        // Compute Atb = Aáµ— * B
-        int incx = 1;
-        int incy = 1;
-        ::dgemv_(&trans, &m, &n, &alpha, A.ptr(), &lda, B.ptr(), &incx, &beta, ws.Atb, &incy);
-
-        // Perform Cholesky factorization of AtA
-        ::dpotrf_(&uplo, &n, ws.AtA, &ldc, &info);
-        if (info != 0) {
-            throw std::runtime_error("Error in DPOTRF: Matrix is not positive definite; info = " + std::to_string(info));
-        }
-
-        // Solve AtA * x = Atb
-        int nrhs = 1;
-        ::dpotrs_(&uplo, &n, &nrhs, ws.AtA, &ldc, ws.Atb, &n, &info);
-        if (info != 0) {
-            throw std::runtime_error("Error in DPOTRS: info = " + std::to_string(info));
-        }
-
-        // Copy solution to weights
-        for (int i = 0; i < n; ++i) {
-            weights[i] = ws.Atb[i];
-        }
-
-    } else {
-        throw std::runtime_error("Unknown algorithm.");
-    }
-}
-
-// Overload without lambda, defaults to lambda = 0.0
-template <typename V1 ,typename V2, typename SVDType>
-void OLS::solve_with_svd(int m_, int n_, V1& B, V2& weights, SVDType& svd) {
-  double rcond = -1.0;
-  double lambda = 0.0;
-  Workspace ws;
-  solve_with_svd_with_workspace(m_, n_, B, weights, svd, ws, rcond, lambda);
-}
-
-// Overload with rcond and lambda
-template <typename V1 ,typename V2, typename SVDType>
-void OLS::solve_with_svd(int m_, int n_, V1& B, V2& weights, SVDType& svd, double rcond, double lambda) {
-  Workspace ws;
-  solve_with_svd_with_workspace(m_, n_, B, weights, svd, ws, rcond, lambda);
-}
-
-// Overload with workspace and default lambda = 0.0
-template <typename V1 ,typename V2, typename SVDType>
-void OLS::solve_with_svd_with_workspace(int m_, int n_, V1& B, V2& weights, SVDType& svd, Workspace& ws) {
-  double rcond = -1.0;
-  double lambda = 0.0;
-  solve_with_svd_with_workspace(m_, n_, B, weights, svd, ws, rcond, lambda);
-}
-
-// Main method with workspace, rcond, and lambda
-template <typename V1 ,typename V2, typename SVDType>
-void OLS::solve_with_svd_with_workspace(int m_, int n_, V1& B, V2& weights, SVDType& svd, Workspace& ws, double rcond, double lambda) {
-  int m = m_;
-  int n = n_;
-
-  int min_mn = std::min(m, n);
-
-  // Check that dimensions of weights match expected size
-  if (weights.size() != static_cast<std::size_t>(n)) {
-    throw std::invalid_argument("weights size is incorrect.");
-  }
-
-  // Check that B has sufficient size
-  int maxmn = std::max(m, n);
-  if (B.size() < static_cast<std::size_t>(maxmn)) {
-    throw std::invalid_argument("B size is insufficient.");
-  }
-
-  // Check if workspace is sufficient; if not, reallocate
-  if (!ws.is_sufficient(m, n, Algorithm::SVD)) {
-    ws.allocate(m, n, Algorithm::SVD);
-  }
-
-  // Copy precomputed SVD components into the workspace
-  ws.U = svd.getU();   // Assume getU() returns a pointer to U (m x min(m,n))
-  ws.S = svd.getS();   // Assume getS() returns a pointer to singular values (min(m,n))
-  ws.VT = svd.getVT(); // Assume getVT() returns a pointer to VT (min(m,n) x n)
-  ws.lambda = lambda;
-
-  // Indicate that we don't own the SVD component pointers
-  ws.owns_U = false;
-  ws.owns_S = false;
-  ws.owns_VT = false;
-
-  // Proceed with the computations
-  if (ws.U_t_b == nullptr) {
-    ws.U_t_b = new double[min_mn];
-  }
-
-  // Compute Uáµ— * B
-  char trans = 'T'; // Transpose U
-  double alpha = 1.0;
-  double beta = 0.0;
-  int incx = 1;
-  int incy = 1;
-  int lda = m;
-  ::dgemv_(&trans, &m, &min_mn, &alpha, ws.U, &lda, B.ptr(), &incx, &beta, ws.U_t_b, &incy);
-
-  // Find the maximum singular value
-  double smax = 0.0;
-  for (int i = 0; i < min_mn; ++i) {
-    if (ws.S[i] > smax) {
-      smax = ws.S[i];
-    }
-  }
-
-  // Threshold for singular values
-  double thresh = rcond > 0.0 ? rcond * smax : 0.0;
+} // namespace models
+} // namespace tadah
 
-  // Apply regularization using 'smax' approach and lambda
-  for (int i = 0; i < min_mn; ++i) {
-    if (ws.S[i] > thresh) {
-      ws.U_t_b[i] = ws.U_t_b[i] * ws.S[i] / (ws.S[i] * ws.S[i] + ws.lambda);
-    } else {
-      ws.U_t_b[i] = 0.0;
-    }
-  }
+#include <tadah/models/ols_impl.h> // Include implementation
 
-  // Compute x = V * (regularized solution)
-  trans = 'T'; // Transpose VT to get V
-  int lda_vt = min_mn;
-  ::dgemv_(&trans, &min_mn, &n, &alpha, ws.VT, &lda_vt, ws.U_t_b, &incx, &beta, weights.ptr(), &incy);
-}
+#endif // TADAH_MODELS_OLS_H
 
-#endif // OLS_H
diff --git a/include/tadah/models/ols_impl.h b/include/tadah/models/ols_impl.h
new file mode 100644
index 0000000..8a0b61a
--- /dev/null
+++ b/include/tadah/models/ols_impl.h
@@ -0,0 +1,311 @@
+#ifndef TADAH_MODELS_OLS_IMPL_H
+#define TADAH_MODELS_OLS_IMPL_H
+
+#include <tadah/core/lapack.h>
+#include <tadah/models/ols.h>
+#include <tadah/models/memory/OLSWorkspace.h>
+#include <tadah/models/memory/SVDWorkspace.h>
+#include <tadah/models/memory/ModelsWorkspaceManager.h>
+#include <stdexcept>
+#include <string>
+#include <algorithm>
+
+namespace tadah {
+namespace models {
+
+// Helper function to convert Algorithm enum to string
+inline std::string to_string_algorithm(Algorithm algo) {
+    switch (algo) {
+        case Algorithm::GELS:
+            return "GELS";
+        case Algorithm::GELSD:
+            return "GELSD";
+        case Algorithm::Pseudoinverse:
+            return "Pseudoinverse";
+        case Algorithm::NormalEquations:
+            return "NormalEquations";
+        default:
+            return "Unknown Algorithm";
+    }
+}
+
+// Overloaded operator<<
+inline std::ostream& operator<<(std::ostream& os, Algorithm algo) {
+    os << to_string_algorithm(algo);
+    return os;
+}
+
+// Constructors and Destructor
+inline OLS::OLS()
+    : workspaceManager_(new memory::ModelsWorkspaceManager()), ownsWorkspaceManager_(true) {}
+
+inline OLS::OLS(memory::IModelsWorkspaceManager& workspaceManager)
+    : workspaceManager_(&workspaceManager), ownsWorkspaceManager_(false) {}
+
+inline OLS::~OLS() {
+    if (ownsWorkspaceManager_ && workspaceManager_) {
+        delete workspaceManager_;
+    }
+}
+
+// Public solve methods
+template <typename M, typename V>
+void OLS::solve(M& A, V& B, V& weights) {
+    double rcond = -1.0;
+    Algorithm algorithm = Algorithm::GELSD; // Default algorithm suitable for all systems
+    solveImpl(A, B, weights, rcond, algorithm);
+}
+
+template <typename M, typename V>
+void OLS::solve(M& A, V& B, V& weights, double rcond, Algorithm algorithm) {
+    solveImpl(A, B, weights, rcond, algorithm);
+}
+
+// Implementation of solve methods
+template <typename M, typename V>
+void OLS::solveImpl(M& A, V& B, V& weights, double rcond, Algorithm algorithm) {
+    int m = A.rows();
+    int n = A.cols();
+    int maxmn = std::max(m, n);
+    int minmn = std::min(m, n);
+
+    bool isUnderdetermined = (m < n);
+
+    // Check if the algorithm is suitable for the system
+    if (isUnderdetermined) {
+        // Underdetermined system
+        if (!(algorithm == Algorithm::GELSD || algorithm == Algorithm::Pseudoinverse)) {
+            throw std::invalid_argument("The chosen algorithm (" + to_string_algorithm(algorithm) +
+                                        ") is not suitable for underdetermined systems (m < n). "
+                                        "Please use GELSD or Pseudoinverse.");
+        }
+    } else {
+        // Overdetermined or square system
+        // All algorithms are applicable, but Cholesky (NormalEquations) may not handle rank-deficient cases well
+        if (algorithm == Algorithm::NormalEquations && m < n) {
+            throw std::invalid_argument("The NormalEquations algorithm is not recommended for systems where m < n.");
+        }
+    }
+
+    // Ensure that B is of sufficient size
+    if (B.size() < static_cast<std::size_t>(m)) {
+        throw std::invalid_argument("B size is insufficient. It must be at least the number of rows of A.");
+    }
+
+    // Get OLS workspace
+    memory::OLSWorkspace* ws = workspaceManager_->getOLSWorkspace(m, n, algorithm);
+
+    // Ensure dimensions match
+    if (weights.size() != static_cast<std::size_t>(n)) {
+        throw std::invalid_argument("weights size is incorrect. It must be of size n (number of columns of A).");
+    }
+
+    // Prepare the right-hand side vector
+    V* b_work_ptr;  // Pointer to the vector to be used as B in computations
+    V b_work;       // Temporary vector if needed
+
+    // If B.size() >= maxmn, we can use B directly
+    if (B.size() >= static_cast<std::size_t>(maxmn)) {
+        b_work_ptr = &B;
+    } else {
+        // Create a temporary vector of size maxmn, initialize with B and zeros
+        b_work.resize(maxmn);
+        for (int i = 0; i < m; ++i) {
+            b_work[i] = B[i];
+        }
+for (int i = m; i < maxmn; ++i) {
+  b_work[i] = 0.0;
+}
+        b_work_ptr = &b_work;
+    }
+
+    int info;
+    if (algorithm == Algorithm::GELSD) {
+        // Use DGELSD for underdetermined or rank-deficient systems
+        int nrhs = 1;
+        double* a = A.ptr();
+        int lda = m;
+        double* b = b_work_ptr->ptr();  // Use the b_work vector
+        int ldb = maxmn;  // Must be at least max(m,n)
+        int rank;
+
+        // Solve using DGELSD
+        ::dgelsd_(&m, &n, &nrhs, a, &lda, b, &ldb, ws->s, &rcond, &rank,
+                  ws->work, &ws->lwork, ws->iwork, &info);
+
+        if (info != 0) {
+            throw std::runtime_error("Error in DGELSD: info = " + std::to_string(info));
+        }
+
+        // Copy solution to weights
+        for (int i = 0; i < n; ++i) {
+            weights[i] = b[i];
+        }
+
+    } else if (algorithm == Algorithm::GELS) {
+        // Use DGELS for overdetermined systems
+        if (isUnderdetermined) {
+            throw std::invalid_argument("GELS algorithm cannot handle underdetermined systems (m < n). Please use GELSD or Pseudoinverse.");
+        }
+        int nrhs = 1;
+        double* a = A.ptr();
+        int lda = m;
+        double* b = b_work_ptr->ptr();
+        int ldb = m;
+        char trans = 'N';
+
+        // Solve using DGELS
+        ::dgels_(&trans, &m, &n, &nrhs, a, &lda, b, &ldb, ws->work, &ws->lwork, &info);
+
+        if (info != 0) {
+            throw std::runtime_error("Error in DGELS: info = " + std::to_string(info));
+        }
+
+        // Copy solution to weights
+        for (int i = 0; i < n; ++i) {
+            weights[i] = b[i];
+        }
+    } else if (algorithm == Algorithm::Pseudoinverse) {
+        // Custom SVD-based pseudoinverse solution
+        // Compute the pseudoinverse using SVD and solve x = A^+ B
+
+        // Copy A into U (since A might be overwritten)
+        /*std::copy(A.ptr(), A.ptr() + m * n, ws->U);*/
+
+        // Perform SVD
+        char jobz = 'S';  // Compute the singular vectors
+        int lda = m;
+        int ldu = m;
+        int ldvt = n;  // For 'S', VT has dimensions (n x min(m, n))
+    ::dgesdd_(&jobz, &m, &n, A.ptr(), &lda, ws->S, ws->U, &ldu,
+              ws->VT, &ldvt, ws->work, &ws->lwork, ws->iwork, &info);
+
+        if (info != 0) {
+            throw std::runtime_error("Error in DGESVD: info = " + std::to_string(info));
+        }
+
+        // Compute U^T B
+        char trans = 'T';
+        double alpha = 1.0;
+        double beta = 0.0;
+        int incx = 1;
+        int incy = 1;
+        ::dgemv_(&trans, &m, &minmn, &alpha, ws->U, &lda, B.ptr(), &incx, &beta, ws->UtB, &incy);
+
+        // Apply threshold to singular values (for handling rank deficiency)
+        double smax = *std::max_element(ws->S, ws->S + minmn);
+        double thresh = (rcond > 0.0) ? rcond * smax : 0.0;
+
+    for (int i = 0; i < minmn; ++i) {
+        if (ws->S[i] > thresh) {
+            ws->UtB[i] /= ws->S[i];
+        } else {
+            ws->UtB[i] = 0.0;
+        }
+    }
+
+        // Compute x = V * (pseudoinverse solution)
+        trans = 'T'; // Use V (not transposed)
+        /*int ldvt = n;*/
+        ::dgemv_(&trans, &minmn, &n, &alpha, ws->VT, &ldvt, ws->UtB, &incx, &beta, weights.ptr(), &incy);
+
+    } else if (algorithm == Algorithm::NormalEquations) {
+        // Solve using normal equations (A^T A x = A^T B)
+        // Compute A^T A
+        char uplo = 'U';
+        char trans = 'T';
+        double alpha = 1.0;
+        double beta = 0.0;
+        int lda = m;
+        int ldc = n;
+
+        ::dsyrk_(&uplo, &trans, &n, &m, &alpha, A.ptr(), &lda, &beta, ws->AtA, &ldc);
+
+        // Compute A^T B
+        int incx = 1;
+        int incy = 1;
+        ::dgemv_(&trans, &m, &n, &alpha, A.ptr(), &lda, B.ptr(), &incx, &beta, ws->Atb, &incy);
+
+        // Perform Cholesky factorization of A^T A
+        ::dpotrf_(&uplo, &n, ws->AtA, &ldc, &info);
+        if (info != 0) {
+            throw std::runtime_error("Error in DPOTRF: Matrix is not positive definite; info = " + std::to_string(info));
+        }
+
+        // Solve A^T A x = A^T B
+        int nrhs = 1;
+        ::dpotrs_(&uplo, &n, &nrhs, ws->AtA, &ldc, ws->Atb, &n, &info);
+        if (info != 0) {
+            throw std::runtime_error("Error in DPOTRS: info = " + std::to_string(info));
+        }
+
+        // Copy solution to weights
+        for (int i = 0; i < n; ++i) {
+            weights[i] = ws->Atb[i];
+        }
+
+    } else {
+        throw std::invalid_argument("Unknown or unsupported algorithm.");
+    }
+
+    // No need to resize B back, since we didn't modify the user's B
+}
+
+// Solve with SVD workspace (if needed)
+template <typename V1, typename V2>
+void OLS::solveWithPseudoinverse(int m, int n, V1& B, V2& weights, memory::SVDWorkspace& svdWorkspace) {
+  double rcond=-1;
+  double lambda=0;
+  OLS::solveWithPseudoinverse(m,n,B,weights,svdWorkspace,rcond,lambda);
+}
+template <typename V1, typename V2>
+void OLS::solveWithPseudoinverse(int m, int n, V1& B, V2& weights, memory::SVDWorkspace& svdWorkspace, double rcond, double lambda) {
+    // Function implementation without default argument
+    memory::OLSWorkspace* ws = workspaceManager_->getOLSWorkspace(m, n, Algorithm::Pseudoinverse);
+
+    // Ensure dimensions match
+    if (weights.size() != static_cast<std::size_t>(n)) {
+        throw std::invalid_argument("weights size is incorrect.");
+    }
+
+    int minmn = std::min(m, n);
+
+    // Use svdWorkspace's U, S, VT directly
+    double* U = svdWorkspace.U;
+    double* S = svdWorkspace.S;
+    double* VT = svdWorkspace.VT;
+
+    // Compute U^T * B
+    char trans = 'T';
+    double alpha = 1.0;
+    double beta = 0.0;
+    int incx = 1;
+    int incy = 1;
+    int lda = m;
+    ::dgemv_(&trans, &m, &minmn, &alpha, U, &lda, B.ptr(), &incx, &beta, ws->UtB, &incy);
+
+    // Find the maximum singular value
+    double smax = *std::max_element(S, S + minmn);
+
+    // Threshold for singular values
+    double thresh = (rcond > 0.0) ? rcond * smax : 0.0;
+
+    // Compute the pseudoinverse solution
+    for (int i = 0; i < minmn; ++i) {
+        if (S[i] > thresh) {
+            ws->UtB[i] *= S[i] / ( S[i] * S[i] + lambda );
+        } else {
+            ws->UtB[i] = 0.0;
+        }
+    }
+
+    // Compute x = V * (pseudoinverse solution)
+    trans = 'T'; // Use V
+    int ldvt = n;
+    ::dgemv_(&trans, &minmn, &n, &alpha, VT, &ldvt, ws->UtB, &incx, &beta, weights.ptr(), &incy);
+}
+
+} // namespace models
+} // namespace tadah
+
+#endif // TADAH_MODELS_OLS_IMPL_H
diff --git a/include/tadah/models/ridge_regression.h b/include/tadah/models/ridge_regression.h
deleted file mode 100644
index 7f208e3..0000000
--- a/include/tadah/models/ridge_regression.h
+++ /dev/null
@@ -1,110 +0,0 @@
-#ifndef RIDGE_REGRESSION_H
-#define RIDGE_REGRESSION_H
-
-#include <tadah/core/core_types.h>
-#include <tadah/core/maths.h>
-#include <tadah/core/lapack.h>
-#include <tadah/models/svd.h>
-#include <algorithm> // For std::max_element
-
-/**
- * @class RidgeRegression
- * @brief Implements ridge regression using Singular Value Decomposition (SVD).
- *
- * Provides methods to perform ridge regression, a type of linear regression that includes
- * regularization to prevent overfitting.
- */
-class RidgeRegression {
-  public:
-    /**
-     * @brief Inverts and multiplies sigma values by lambda for regularization,
-     *        considering the rcond threshold.
-     *
-     * Performs an element-wise operation on the sigma values to handle regularization.
-     * Singular values below the threshold are discarded (set to zero).
-     *
-     * @tparam V Vector type for sigma and result.
-     * @param sigma Input vector of sigma values.
-     * @param result Output vector for inverted and multiplied results.
-     * @param n Size of the vectors.
-     * @param lambda Regularization parameter.
-     * @param rcond Reciprocal of the condition number; threshold for singular values.
-     */
-    template <typename V>
-    static void invertMultiplySigmaLambdaRcond(const V& sigma, V& result, int n, double lambda, double rcond) {
-      // Find the maximum singular value
-      double smax = 0.0;
-      for (int i = 0; i < n; ++i) {
-        if (sigma[i] > smax) {
-          smax = sigma[i];
-        }
-      }
-
-      // Threshold for singular values
-      double thresh = rcond > 0.0 ? rcond * smax : 0.0;
-
-      for (int i = 0; i < n; ++i) {
-        if (sigma[i] > thresh) {
-          result[i] = sigma[i] / (sigma[i] * sigma[i] + lambda);
-        } else {
-          // Singular value is too small; set result to zero
-          result[i] = 0.0;
-        }
-      }
-    }
-
-    /**
-     * @brief Solves the ridge regression problem using SVD components,
-     *        incorporating rcond for singular value thresholding.
-     *
-     * Computes the weights that minimize the regularized least squares error.
-     *
-     * @tparam V Vector type for b.
-     * @tparam W Vector type for weights.
-     * @param svd SVD object containing decomposed matrix components.
-     * @param b Vector of target values.
-     * @param weights Output vector for computed weights.
-     * @param lambda Regularization parameter.
-     * @param rcond Reciprocal of the condition number; threshold for singular values.
-     */
-    template <typename V, typename W>
-    static void solve(const SVD &svd, V b, W &weights, double lambda, double rcond) {
-      double *U = svd.getU();    // Matrix U from SVD (m x n)
-      double *s = svd.getS();    // Singular values (as a vector)
-      double *VT = svd.getVT();  // Matrix V^T from SVD (n x n)
-
-      int m = svd.shapeA().first;
-      int n = svd.shapeA().second;
-
-      // Allocate memory for intermediate calculations
-      double *d = new double[n];   // For inverted singular values
-      double *UTb = new double[n]; // For U^T * b
-
-      // Step 1: Compute U^T * b
-      char trans = 'T';
-      double alpha_ = 1.0;  // Scalar for multiplication
-      double beta_ = 0.0;   // Scalar for addition
-      int incx = 1;
-      int incy = 1;
-      dgemv_(&trans, &m, &n, &alpha_, U, &m, b.ptr(), &incx, &beta_, UTb, &incy);
-
-      // Step 2: Invert and multiply sigma values, applying rcond threshold
-      invertMultiplySigmaLambdaRcond(s, d, n, lambda, rcond);
-
-      // Element-wise multiplication: d[i] = d[i] * (U^T * b)[i]
-      for (int i = 0; i < n; ++i) {
-        d[i] *= UTb[i];
-      }
-
-      // Step 3: Compute weights = V * d
-      trans = 'T'; // Since VT is V^T, transposing VT gives V
-      weights.resize(n);
-      dgemv_(&trans, &n, &n, &alpha_, VT, &n, d, &incx, &beta_, weights.ptr(), &incy);
-
-      // Cleanup dynamic memory
-      delete[] d;
-      delete[] UTb;
-    }
-};
-
-#endif // RIDGE_REGRESSION_H
diff --git a/include/tadah/models/svd.h b/include/tadah/models/svd.h
index 707c20d..45d3f4b 100644
--- a/include/tadah/models/svd.h
+++ b/include/tadah/models/svd.h
@@ -1,285 +1,127 @@
-#ifndef SVD_H
-#define SVD_H
+#ifndef TADAH_MODELS_SVD_H
+#define TADAH_MODELS_SVD_H
 
-#include <tadah/core/core_types.h>
-#include <tadah/core/maths.h>
-#include <tadah/core/lapack.h>
+#include <tadah/models/memory/IModelsWorkspaceManager.h>
+#include <tadah/models/memory/SVDWorkspace.h>
+
+namespace tadah {
+namespace models {
 
 /**
  * @class SVD
- * @brief Computes the Singular Value Decomposition (SVD) of a matrix.
+ * @brief Provides functionality for computing the Singular Value Decomposition (SVD) of a matrix.
+ *
+ * The SVD class computes the SVD of a given matrix \f$ A \f$, decomposing it into
+ * \f$ A = U \Sigma V^\top \f$, where \f$ U \f$ and \f$ V \f$ are orthogonal matrices,
+ * and \f$ \Sigma \f$ is a diagonal matrix containing the singular values.
  *
- * Decomposes a matrix into U, S, and V^T such that M = USV^T. Provides methods for accessing
- * components and predicting the effect of this decomposition.
+ * **Usage:**
+ * ```cpp
+ * tadah::models::SVD svd;
+ * svd.compute(A);
+ * double* U = svd.getU();
+ * double* S = svd.getS();
+ * double* VT = svd.getVT();
+ * ```
  */
 class SVD {
-  private:
-    char jobz = 'S'; ///< Compute the first min(M,N) columns of U
-    int m, n, lda, ldu, ldvt, ucol, lwork, info;
-
-    // Arrays to store SVD components
-    double *a = nullptr;
-    double *s = nullptr;
-    double *u = nullptr;
-    double *vt = nullptr;
-    double *work = nullptr;
-    int *iwork = nullptr;
-
-    bool work_alloc = false;
-
-  public:
-    /**
-     * @brief Destructor to free allocated resources.
-     */
-    ~SVD() {
-      delete[] s;
-      delete[] u;
-      delete[] vt;
-
-      if (work_alloc)
-        delete[] iwork;
-
-      if (work_alloc && lwork > 0)
-        delete[] work;
-    }
-
-    /**
-     * @brief Copy constructor for SVD class.
-     * 
-     * Performs a deep copy of another SVD instance.
-     * 
-     * @param other The SVD instance to copy.
-     */
-    SVD(const SVD &other) {
-      jobz = other.jobz;
-      m = other.m;
-      n = other.n;
-      lda = other.lda;
-      ldu = other.ldu;
-      ldvt = other.ldvt;
-      ucol = other.ucol;
-      lwork = other.lwork;
-      info = other.info;
-
-      // Allocate memory
-      s = new double[sizeS()];
-      u = new double[sizeU()];
-      vt = new double[sizeVT()];
-
-      // No need for these arrays in a copy
-      a = nullptr;
-      work = nullptr;
-      iwork = nullptr;
-
-      work_alloc = false;
-
-      // Hard copy the data
-      std::memcpy(s, other.s, sizeS() * sizeof(double));
-      std::memcpy(u, other.u, sizeU() * sizeof(double));
-      std::memcpy(vt, other.vt, sizeVT() * sizeof(double));
-    }
-
+public:
     /**
-     * @brief Constructs the SVD of the given matrix.
+     * @brief Default constructor.
      *
-     * @param Phi Input matrix to decompose.
-     */
-    SVD(Matrix &Phi)
-      : m(Phi.rows()),
-        n(Phi.cols()),
-        lda(Phi.rows()),
-        ldu(Phi.rows()),
-        ldvt(std::min(m, n)),
-        ucol(std::min(m, n)),
-        a(&Phi.data()[0]),
-        s(new double[sizeS()]),
-        u(new double[sizeU()]),
-        vt(new double[sizeVT()]),
-        iwork(new int[8 * std::min(m, n)]),
-        work_alloc(true)
-    {
-      lwork = -1; // Query optimal work size
-      double wkopt;
-      dgesdd_(&jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
-          &wkopt, &lwork, iwork, &info);
-      lwork = (int)wkopt;
-      work = new double[lwork];
-      dgesdd_(&jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt,
-          work, &lwork, iwork, &info);
-
-    //   // Find the maximum singular value
-    //   double smax = 0.0;
-    //   for (int i = 0; i < n; ++i) {
-    //     if (s[i] > smax) {
-    //       smax = s[i];
-    //     }
-    //   }
-    //
-    //   // Threshold for singular values
-    //   double thresh = rcond > 0.0 ? rcond * smax : 0.0;
-    //
-    // for (int i = 0; i < n; ++i) {
-    //   if (s[i] < thresh) {
-    //     // Singular value is too small; set result to zero
-    //     s[i] = 0.0;
-    //   }
-    // }
-    }
-
-    /**
-     * @brief Returns the computational info from the LAPACK routine.
-     * @return Info status.
-     */
-    int get_info() const {
-      return info;
-    }
-
-    /**
-     * @brief Accessor for the original matrix A.
-     * @return Pointer to matrix A.
+     * Creates an instance of SVD with its own workspace manager.
      */
-    double *getA() const {
-      return a;
-    }
+    SVD();
 
     /**
-     * @brief Accessor for the U matrix.
-     * @return Pointer to matrix U.
-     */
-    double *getU() const {
-      return u;
-    }
-
-    /**
-     * @brief Accessor for the V^T matrix.
-     * @return Pointer to matrix V^T.
-     */
-    double *getVT() const {
-      return vt;
-    }
-
-    /**
-     * @brief Accessor for the singular values matrix S.
-     * @return Pointer to singular values.
-     */
-    double *getS() const {
-      return s;
-    }
-
-    /**
-     * @brief Returns the size of matrix A.
-     * @return Size of matrix A.
+     * @brief Constructor with external workspace manager.
+     *
+     * @param workspaceManager Reference to an external IModelsWorkspaceManager.
      */
-    size_t sizeA() const {
-      return m * n;
-    }
+    explicit SVD(memory::IModelsWorkspaceManager& workspaceManager);
 
     /**
-     * @brief Returns the size of matrix U.
-     * @return Size of matrix U.
+     * @brief Destructor.
+     *
+     * Deletes the workspace manager if owned.
      */
-    size_t sizeU() const {
-      return ldu * ucol;
-    }
+    ~SVD();
 
     /**
-     * @brief Returns the size of matrix V^T.
-     * @return Size of matrix V^T.
+     * @brief Compute the SVD of a given matrix A.
+     *
+     * @tparam M Type of the matrix A.
+     * @param A The matrix A to decompose.
+     *
+     * **Example:**
+     * ```cpp
+     * tadah::core::Matrix A; // Initialize your matrix
+     * tadah::models::SVD svd;
+     * svd.compute(A);
+     * ```
      */
-    size_t sizeVT() const {
-      return ldvt * n;
-    }
+    template <typename M>
+    void compute(M& A);
 
     /**
-     * @brief Returns the number of singular values.
-     * @return Number of singular values.
+     * @brief Get the U matrix from the SVD.
+     *
+     * @return Pointer to the U matrix (size m x min(m, n)).
      */
-    size_t sizeS() const {
-      return std::min(m, n);
-    }
+    double* getU() const;
 
     /**
-     * @brief Returns the dimensions of matrix A.
-     * @return Pair of dimensions (rows, cols).
+     * @brief Get the singular values from the SVD.
+     *
+     * @return Pointer to the singular values array (size min(m, n)).
      */
-    std::pair<size_t, size_t> shapeA() const {
-      return std::pair<size_t, size_t>(m, n);
-    }
+    double* getS() const;
 
     /**
-     * @brief Returns the dimensions of matrix U.
-     * @return Pair of dimensions (rows, cols).
+     * @brief Get the V^T matrix from the SVD.
+     *
+     * @return Pointer to the V^T matrix (size min(m, n) x n).
      */
-    std::pair<size_t, size_t> shapeU() const {
-      return std::pair<size_t, size_t>(ldu, ucol);
-    }
+    double* getVT() const;
 
     /**
-     * @brief Returns the dimensions of matrix V^T.
-     * @return Pair of dimensions (rows, cols).
+     * @brief Get the SVDWorkspace used in the computation.
+     *
+     * @return Reference to the SVDWorkspace.
      */
-    std::pair<size_t, size_t> shapeVT() const {
-      return std::pair<size_t, size_t>(ldvt, n);
-    }
+    memory::SVDWorkspace& getSVDWorkspace();
 
     /**
-     * @brief Returns the dimensions of the singular values.
-     * @return Pair of dimensions (count, 1).
+     * @brief Get the number of rows of the original matrix.
+     *
+     * @return Number of rows.
      */
-    std::pair<size_t, size_t> shapeS() const {
-      return std::pair<size_t, size_t>(std::min(m, n), 1);
-    }
+    int rows() const;
 
     /**
-     * @brief Predicts the result of multiplying U, S, V^T with a vector.
+     * @brief Get the number of columns of the original matrix.
      *
-     * @tparam T Type of vector w.
-     * @param w Vector to be transformed.
-     * @return Transformed vector.
+     * @return Number of columns.
      */
-    template <typename T>
-    T predict(const T &w) const {
-      double *t = new double[shapeA().first];
-      double *x = new double[shapeA().second];
+    int cols() const;
 
-      T res = predict(w, t, x);
+private:
+    memory::IModelsWorkspaceManager* workspaceManager_; ///< Pointer to the workspace manager.
+    bool ownsWorkspaceManager_; ///< Indicates if this instance owns the workspace manager.
+    memory::SVDWorkspace* svdWorkspace_; ///< Pointer to the SVDWorkspace.
 
-      delete[] t;
-      delete[] x;
-      return res;
-    }
+    int m_; ///< Number of rows in the original matrix.
+    int n_; ///< Number of columns in the original matrix.
 
-    /**
-     * @brief Predicts the result with intermediate storage.
-     *
-     * @tparam T Type of vector w.
-     * @param w Vector to be transformed.
-     * @param t Temporary storage for transformation.
-     * @param x Temporary storage for transformation.
-     * @return Transformed vector.
-     */
-    template <typename T>
-    T predict(const T &w, double *t, double *x) const {
-      const double *w_ptr = &w.data()[0];
-      char trans = 'N';
-      int m = static_cast<int>(shapeA().first);
-      int n = static_cast<int>(shapeA().second);
-      double alpha_ = 1.0;
-      double beta_ = 0.0;
-      int incx = 1;
+    // Disable copying
+    SVD(const SVD&) = delete;
+    SVD& operator=(const SVD&) = delete;
+};
 
-      // 1. V^T w
-      dgemv_(&trans, &n, &n, &alpha_, vt, &n, w_ptr, &incx, &beta_, x, &incx);
+} // namespace models
+} // namespace tadah
 
-      // 2. S (V^T w)
-      for (int i = 0; i < n; ++i) {
-        x[i] *= s[i];
-      }
+#include <tadah/models/svd_impl.h> // Include implementation
 
-      // 3. U (S V^T w)
-      dgemv_(&trans, &m, &n, &alpha_, u, &m, x, &incx, &beta_, t, &incx);
-      return T(t, m);
-    }
-};
+#endif // TADAH_MODELS_SVD_H
 
-#endif // SVD_H
diff --git a/include/tadah/models/svd_impl.h b/include/tadah/models/svd_impl.h
new file mode 100644
index 0000000..a7fb017
--- /dev/null
+++ b/include/tadah/models/svd_impl.h
@@ -0,0 +1,79 @@
+#ifndef TADAH_MODELS_SVD_IMPL_H
+#define TADAH_MODELS_SVD_IMPL_H
+
+#include <tadah/models/svd.h>
+#include <tadah/models/memory/ModelsWorkspaceManager.h>
+#include <tadah/core/lapack.h>
+#include <algorithm>
+#include <stdexcept>
+
+namespace tadah {
+namespace models {
+
+inline SVD::SVD()
+    : workspaceManager_(new memory::ModelsWorkspaceManager()), ownsWorkspaceManager_(true), svdWorkspace_(nullptr), m_(0), n_(0) {}
+
+inline SVD::SVD(memory::IModelsWorkspaceManager& workspaceManager)
+    : workspaceManager_(&workspaceManager), ownsWorkspaceManager_(false), svdWorkspace_(nullptr), m_(0), n_(0) {}
+
+inline SVD::~SVD() {
+    if (ownsWorkspaceManager_ && workspaceManager_) {
+        delete workspaceManager_;
+    }
+}
+
+template <typename M>
+void SVD::compute(M& A) {
+    int m = A.rows();
+    int n = A.cols();
+    m_ = m; // Store the number of rows
+    n_ = n; // Store the number of columns
+    svdWorkspace_ = workspaceManager_->getSVDWorkspace(m, n);
+
+    int minMN = std::min(m, n);
+
+    // Prepare parameters for LAPACK dgesdd_
+    char jobz = 'S'; // Compute singular values and vectors
+    int lda = m;
+    int ldu = m;
+    int ldvt = minMN;
+    int info;
+
+    // Perform SVD
+    ::dgesdd_(&jobz, &m, &n, A.ptr(), &lda, svdWorkspace_->S, svdWorkspace_->U, &ldu,
+              svdWorkspace_->VT, &ldvt, svdWorkspace_->work, &svdWorkspace_->lwork, svdWorkspace_->iwork, &info);
+
+    if (info != 0) {
+        throw std::runtime_error("Error in DGESDD during SVD computation: info = " + std::to_string(info));
+    }
+}
+
+inline double* SVD::getU() const {
+    return svdWorkspace_->U;
+}
+
+inline double* SVD::getS() const {
+    return svdWorkspace_->S;
+}
+
+inline double* SVD::getVT() const {
+    return svdWorkspace_->VT;
+}
+
+inline memory::SVDWorkspace& SVD::getSVDWorkspace() {
+    return *svdWorkspace_;
+}
+
+inline int SVD::rows() const {
+    return m_;
+}
+
+inline int SVD::cols() const {
+    return n_;
+}
+
+} // namespace models
+} // namespace tadah
+
+#endif // TADAH_MODELS_SVD_IMPL_H
+
diff --git a/src/ModelsWorkspaceManager.cpp b/src/ModelsWorkspaceManager.cpp
new file mode 100644
index 0000000..485c1d1
--- /dev/null
+++ b/src/ModelsWorkspaceManager.cpp
@@ -0,0 +1,54 @@
+#include <tadah/models/memory/ModelsWorkspaceManager.h>
+#include <tadah/models/memory/OLSWorkspace.h>
+#include <tadah/models/memory/SVDWorkspace.h>
+
+namespace tadah {
+namespace models {
+namespace memory {
+
+ModelsWorkspaceManager::ModelsWorkspaceManager()
+    : olsWorkspace_(nullptr), svdWorkspace_(nullptr) {}
+
+ModelsWorkspaceManager::~ModelsWorkspaceManager() {
+    releaseOLSWorkspace(olsWorkspace_);
+    releaseSVDWorkspace(svdWorkspace_);
+}
+
+OLSWorkspace* ModelsWorkspaceManager::getOLSWorkspace(int m, int n, tadah::models::Algorithm algorithm) {
+    if (!olsWorkspace_ || !olsWorkspace_->isSufficient(m, n, algorithm)) {
+        releaseOLSWorkspace(olsWorkspace_);
+        olsWorkspace_ = new OLSWorkspace();
+        olsWorkspace_->allocate(m, n, algorithm);
+    }
+    return olsWorkspace_;
+}
+
+void ModelsWorkspaceManager::releaseOLSWorkspace(OLSWorkspace* workspace) {
+    if (workspace) {
+        delete workspace;
+        workspace = nullptr;
+    }
+    olsWorkspace_ = nullptr;
+}
+
+SVDWorkspace* ModelsWorkspaceManager::getSVDWorkspace(int m, int n) {
+    if (!svdWorkspace_ || !svdWorkspace_->isSufficient(m, n)) {
+        releaseSVDWorkspace(svdWorkspace_);
+        svdWorkspace_ = new SVDWorkspace();
+        svdWorkspace_->allocate(m, n);
+    }
+    return svdWorkspace_;
+}
+
+void ModelsWorkspaceManager::releaseSVDWorkspace(SVDWorkspace* workspace) {
+    if (workspace) {
+        delete workspace;
+        workspace = nullptr;
+    }
+    svdWorkspace_ = nullptr;
+}
+
+} // namespace memory
+} // namespace models
+} // namespace tadah
+
diff --git a/src/OLSWorkspace.cpp b/src/OLSWorkspace.cpp
new file mode 100644
index 0000000..aab90d4
--- /dev/null
+++ b/src/OLSWorkspace.cpp
@@ -0,0 +1,167 @@
+#include <algorithm>
+#include <cmath>
+#include <stdexcept>
+#include <tadah/core/lapack.h>
+#include <tadah/models/memory/OLSWorkspace.h>
+
+namespace tadah {
+namespace models {
+namespace memory {
+
+OLSWorkspace::OLSWorkspace()
+    : work(nullptr), lwork(-1), m(0), n(0),
+      algo(tadah::models::Algorithm::GELSD), s(nullptr), iwork(nullptr),
+      lambda(0.0), U(nullptr), S(nullptr), VT(nullptr), iworkSVD(nullptr),
+      UtB(nullptr), AtA(nullptr), Atb(nullptr) {}
+
+OLSWorkspace::~OLSWorkspace() {
+  // Deallocate all allocated memory
+  delete[] U;
+  delete[] S;
+  delete[] VT;
+  delete[] s;
+  delete[] iwork;
+  delete[] UtB;
+  delete[] iworkSVD;
+  delete[] work;
+  delete[] AtA;
+  delete[] Atb;
+
+  // Reset pointers
+  U = nullptr;
+  S = nullptr;
+  VT = nullptr;
+  s = nullptr;
+  iwork = nullptr;
+  UtB = nullptr;
+  iworkSVD = nullptr;
+  work = nullptr;
+  AtA = nullptr;
+  Atb = nullptr;
+}
+
+void OLSWorkspace::allocate(int m_, int n_,
+                            tadah::models::Algorithm algorithm) {
+  // Release existing memory
+  delete[] U;
+  delete[] S;
+  delete[] VT;
+  delete[] s;
+  delete[] iwork;
+  delete[] UtB;
+  delete[] iworkSVD;
+  delete[] work;
+  delete[] AtA;
+  delete[] Atb;
+
+  // Reset pointers
+  U = nullptr;
+  S = nullptr;
+  VT = nullptr;
+  s = nullptr;
+  iwork = nullptr;
+  UtB = nullptr;
+  iworkSVD = nullptr;
+  work = nullptr;
+  AtA = nullptr;
+  Atb = nullptr;
+
+  m = m_;
+  n = n_;
+  algo = algorithm;
+
+  int info = 0;
+
+  if (algorithm == tadah::models::Algorithm::GELSD) {
+    // Allocate workspace for DGELSD
+    int minmn = std::min(m, n);
+    int maxmn = std::max(m, n);
+    int nlvl = std::max(0, static_cast<int>(std::log2(minmn / 25.0)) + 1);
+
+    // Allocate s (singular values)
+    s = new double[minmn];
+
+    // Allocate iwork
+    int liwork = 3 * minmn * nlvl + 11 * minmn;
+    iwork = new int[liwork];
+
+    // Workspace query
+    double workQuery;
+    int lworkQuery = -1;
+    int nrhs = 1;
+    double *aDummy = nullptr;
+    double *bDummy = nullptr;
+    double rcondDummy = -1.0;
+    int rankDummy;
+    ::dgelsd_(&m, &n, &nrhs, aDummy, &m, bDummy, &maxmn, s, &rcondDummy,
+              &rankDummy, &workQuery, &lworkQuery, iwork, &info);
+    lwork = static_cast<int>(workQuery);
+    work = new double[lwork];
+  } else if (algorithm == tadah::models::Algorithm::GELS) {
+    // Allocate workspace for DGELS
+    double workQuery;
+    int lworkQuery = -1;
+    int nrhs = 1;
+    double *aDummy = nullptr;
+    double *bDummy = nullptr;
+    char trans = 'N';
+    ::dgels_(&trans, &m, &n, &nrhs, aDummy, &m, bDummy, &m, &workQuery,
+             &lworkQuery, &info);
+    lwork = static_cast<int>(workQuery);
+    work = new double[lwork];
+  } else if (algorithm == tadah::models::Algorithm::Pseudoinverse) {
+    // Allocate workspace for internal SVD computation in OLS
+
+    int minMN = std::min(m, n);
+
+    // Allocate U, S, VT
+    U = new double[m * minMN];  // U is m x min(m,n)
+    S = new double[minMN];      // Singular values
+    VT = new double[minMN * n]; // VT is min(m,n) x n
+
+    // Workspace for dgesdd_
+    iworkSVD = new int[8 * minMN];
+
+    // Query for optimal workspace size
+    char jobz = 'S';
+    int lda = m;
+    int ldu = m;
+    int ldvt = minMN;
+    lwork = -1;
+    double wkopt;
+    double *aDummy = nullptr;
+    int infoQuery;
+    ::dgesdd_(&jobz, &m, &n, aDummy, &lda, S, U, &ldu, VT, &ldvt, &wkopt,
+              &lwork, iworkSVD, &infoQuery);
+
+    if (infoQuery != 0) {
+      throw std::runtime_error("Error in DGESDD (workspace query): info = " +
+                               std::to_string(infoQuery));
+    }
+
+    lwork = static_cast<int>(wkopt);
+    work = new double[lwork];
+
+    // Allocate UtB
+    UtB = new double[minMN];
+  } else if (algorithm == tadah::models::Algorithm::NormalEquations) {
+    // Allocate space for AtA (n x n) and Atb (n)
+    AtA = new double[n * n];
+    Atb = new double[n];
+    // No additional workspace required
+    work = nullptr;
+    lwork = 0;
+  } else {
+    throw std::invalid_argument(
+        "Invalid algorithm specified in OLSWorkspace::allocate().");
+  }
+}
+
+bool OLSWorkspace::isSufficient(int m_, int n_,
+                                tadah::models::Algorithm algorithm) const {
+  return (m == m_ && n == n_ && algo == algorithm);
+}
+
+} // namespace memory
+} // namespace models
+} // namespace tadah
diff --git a/src/SVDWorkspace.cpp b/src/SVDWorkspace.cpp
new file mode 100644
index 0000000..172eff4
--- /dev/null
+++ b/src/SVDWorkspace.cpp
@@ -0,0 +1,76 @@
+#include <tadah/models/memory/SVDWorkspace.h>
+#include <tadah/core/lapack.h>
+#include <algorithm>
+#include <cmath>
+#include <stdexcept>
+
+namespace tadah {
+namespace models {
+namespace memory {
+
+SVDWorkspace::SVDWorkspace()
+    : m(0), n(0),
+      U(nullptr), S(nullptr), VT(nullptr),
+      work(nullptr), lwork(-1), iwork(nullptr) {}
+
+SVDWorkspace::~SVDWorkspace() {
+    delete[] U;
+    delete[] S;
+    delete[] VT;
+    delete[] iwork;
+    delete[] work;
+    U = nullptr; S = nullptr; VT = nullptr;
+    iwork = nullptr; work = nullptr;
+}
+
+void SVDWorkspace::allocate(int m_, int n_) {
+    // Release existing memory
+    delete[] U;
+    delete[] S;
+    delete[] VT;
+    delete[] iwork;
+    delete[] work;
+
+    U = nullptr; S = nullptr; VT = nullptr;
+    iwork = nullptr; work = nullptr;
+
+    m = m_;
+    n = n_;
+
+    int minMN = std::min(m, n);
+
+    // Allocate U, S, VT
+    U = new double[m * minMN];      // U is m x min(m,n)
+    S = new double[minMN];          // Singular values
+    VT = new double[minMN * n];     // VT is min(m,n) x n
+
+    // Workspace for dgesdd_
+    iwork = new int[8 * minMN];
+
+    // Query for optimal workspace size
+    char jobz = 'S';
+    int lda = m;
+    int ldu = m;
+    int ldvt = minMN;
+    lwork = -1;
+    double wkopt;
+    double* aDummy = nullptr;
+    int infoQuery;
+    ::dgesdd_(&jobz, &m, &n, aDummy, &lda, S, U, &ldu, VT, &ldvt, &wkopt, &lwork, iwork, &infoQuery);
+
+    if (infoQuery != 0) {
+        throw std::runtime_error("Error in DGESDD (workspace query in SVDWorkspace): info = " + std::to_string(infoQuery));
+    }
+
+    lwork = static_cast<int>(wkopt);
+    work = new double[lwork];
+}
+
+bool SVDWorkspace::isSufficient(int m_, int n_) const {
+    return (m == m_ && n == n_);
+}
+
+} // namespace memory
+} // namespace models
+} // namespace tadah
+
diff --git a/src/ols.cpp b/src/ols.cpp
index 52bc49f..e69de29 100644
--- a/src/ols.cpp
+++ b/src/ols.cpp
@@ -1,165 +0,0 @@
-#include <tadah/models/ols.h>
-#include <cmath>
-#include <stdexcept>
-#include <cstdlib>    // For malloc and free
-#include <algorithm>  // For std::max
-
-// Implementation of non-template methods
-
-// Constructor
-OLS::Workspace::Workspace()
-  : work(nullptr), lwork(-1), m(0), n(0), algo(Algorithm::GELSD),
-  s(nullptr), iwork(nullptr),
-  lambda(0.0), U(nullptr), S(nullptr), VT(nullptr), iwork_svd(nullptr), U_t_b(nullptr),
-  owns_U(false), owns_S(false), owns_VT(false),
-  AtA(nullptr), Atb(nullptr) {
-}
-
-// Destructor
-OLS::Workspace::~Workspace() {
-  release();
-}
-
-// allocate method
-void OLS::Workspace::allocate(int m_, int n_, Algorithm algorithm) {
-  release(); // Release any existing memory
-
-  m = m_;
-  n = n_;
-  algo = algorithm;
-  lambda = 0.0; // Reset lambda
-
-  int info;
-
-  if (algorithm == Algorithm::GELSD) {
-    // Allocate workspace for DGELSD
-    int minmn = std::min(m, n);
-    int maxmn = std::max(m, n);
-    int nlvl = std::max(0, static_cast<int>(std::log2(minmn / 25.0)) + 1);
-
-    // Allocate s (singular values)
-    s = new double[minmn];
-
-    // Allocate iwork
-    int liwork = 3 * minmn * nlvl + 11 * minmn;
-    iwork = new int[liwork];
-
-    // Workspace query
-    double work_query;
-    int lwork_query = -1;
-    int nrhs = 1;
-    double* a_dummy = nullptr;
-    double* b_dummy = nullptr;
-    ::dgelsd_(&m, &n, &nrhs, a_dummy, &m, b_dummy, &maxmn, s, nullptr,
-              nullptr, &work_query, &lwork_query, iwork, &info);
-    lwork = static_cast<int>(work_query);
-    work = new double[lwork];
-  } else if (algorithm == Algorithm::GELS) {
-    // Allocate workspace for DGELS
-    double work_query;
-    int lwork_query = -1;
-    int nrhs = 1;
-    double* a_dummy = nullptr;
-    double* b_dummy = nullptr;
-    char trans = 'N';
-    ::dgels_(&trans, &m, &n, &nrhs, a_dummy, &m, b_dummy, &m, &work_query, &lwork_query, &info);
-    lwork = static_cast<int>(work_query);
-    work = new double[lwork];
-  } else if (algorithm == Algorithm::SVD) {
-
-    // Allocate U, S, VT, and iwork_svd
-    int min_mn = std::min(m, n);
-
-    // Allocate U, S, VT regardless of ownership flags (they should be owned in this case)
-    if (U == nullptr) {
-      U = new double[m * min_mn];      // U is m x min(m,n)
-      owns_U = true;
-    }
-    if (S == nullptr) {
-      S = new double[min_mn];          // Singular values
-      owns_S = true;
-    }
-    if (VT == nullptr) {
-      VT = new double[min_mn * n];     // VT is min(m,n) x n
-      owns_VT = true;
-    }
-    iwork_svd = new int[8 * min_mn];     // Workspace for dgesdd_
-
-    // Query for optimal workspace size
-    char jobz = 'S';
-    int lda = m;
-    int ldu = m;
-    int ldvt = min_mn;
-    lwork = -1;
-    double wkopt;
-    int info_query;
-    double* a_dummy = nullptr;
-    ::dgesdd_(&jobz, &m, &n, a_dummy, &lda, S, U, &ldu, VT, &ldvt, &wkopt, &lwork, iwork_svd, &info_query);
-
-    if (info_query != 0) {
-      throw std::runtime_error("Error in DGESDD (workspace query): info = " + std::to_string(info_query));
-    }
-
-    lwork = static_cast<int>(wkopt);
-    work = new double[lwork];
-
-  } else if (algorithm == Algorithm::Cholesky) {
-    // Allocate space for AtA (n x n) and Atb (n)
-    AtA = new double[n * n];
-    Atb = new double[n];
-
-    // No additional workspace required for Cholesky decomposition
-    work = nullptr;
-    lwork = 0;
-  } else {
-    throw std::invalid_argument("Invalid algorithm specified in Workspace::allocate().");
-  }
-}
-
-// release method
-void OLS::Workspace::release() {
-  if (owns_U && U != nullptr) {
-    delete[] U;
-    U = nullptr;
-    owns_U = false;
-  }
-  if (owns_S && S != nullptr) {
-    delete[] S;
-    S = nullptr;
-    owns_S = false;
-  }
-  if (owns_VT && VT != nullptr) {
-    delete[] VT;
-    VT = nullptr;
-    owns_VT = false;
-  }
-  delete[] s;
-  delete[] iwork;
-  delete[] U_t_b;
-  delete[] iwork_svd;
-  delete[] work;
-  delete[] AtA;
-  delete[] Atb;
-
-  s = nullptr;
-  iwork = nullptr;
-  U_t_b = nullptr;
-  iwork_svd = nullptr;
-  work = nullptr;
-  AtA = nullptr;
-  Atb = nullptr;
-
-  lwork = -1;
-  m = 0;
-  n = 0;
-  algo = Algorithm::GELSD;
-  lambda = 0.0;
-}
-
-// is_sufficient method
-bool OLS::Workspace::is_sufficient(int m_, int n_, Algorithm algorithm) const {
-  if (m != m_ || n != n_ || algo != algorithm) {
-    return false;
-  }
-  return true;
-}
diff --git a/tests/test_ols.cpp b/tests/test_ols.cpp
index 98a34d0..4e0f241 100644
--- a/tests/test_ols.cpp
+++ b/tests/test_ols.cpp
@@ -1,368 +1,354 @@
+// test_ols.cpp
+
 #include "catch2/catch.hpp"
-#include <tadah/core/lapack.h>
-#include <tadah/core/maths.h>
 #include <tadah/models/ols.h>
-#include <iomanip>
+#include <tadah/models/svd.h>
+#include <tadah/models/memory/ModelsWorkspaceManager.h>
+#include <tadah/core/maths.h>
+#include <cmath>
+#include <algorithm>
+#include <stdexcept>
+#include <vector>
 
 TEST_CASE("Testing OLS") {
-  using Matrix = Matrix_T<>; // Replace with your actual matrix type.
-  using Vector = aed_type;   // Replace with your actual vector type.
-
-  // Helper function to compare vectors
-  auto vectors_are_close = [](const Vector& v1, const Vector& v2, double tolerance) {
-    REQUIRE(v1.size() == v2.size());
-    for (size_t i = 0; i < v1.size(); ++i) {
-      REQUIRE_THAT(v1[i], Catch::Matchers::WithinAbs(v2[i], tolerance));
+    using Matrix = Matrix_T<>;
+    using Vector = aed_type;
+
+    // Helper function to compare vectors within a given tolerance
+    auto vectorsAreClose = [](const Vector& v1, const Vector& v2, double tolerance) {
+        REQUIRE(v1.size() == v2.size());
+        for (size_t i = 0; i < v1.size(); ++i) {
+            REQUIRE(std::abs(v1[i] - v2[i]) < tolerance);
+        }
+    };
+
+    // Prepare test data
+    SECTION("Solve with square matrix") {
+        // Initialize A (column-major order)
+        Matrix A(3, 3);
+        A(0,0) = 1; A(1,0) = 4; A(2,0) = 7;   // First column
+        A(0,1) = 2; A(1,1) = 5; A(2,1) = 8;   // Second column
+        A(0,2) = 3; A(1,2) = 6; A(2,2) = 10;  // Third column
+
+        Vector expectedWeights = {1, 2, 3};
+
+        // Compute B = A * expectedWeights
+        Vector B = A * expectedWeights;
+
+        // Create an instance of OLS
+        tadah::models::OLS ols;
+
+        // Test different algorithms
+        std::vector<tadah::models::Algorithm> algorithms = {
+            tadah::models::Algorithm::GELSD,
+            tadah::models::Algorithm::GELS,
+            tadah::models::Algorithm::NormalEquations,
+            tadah::models::Algorithm::Pseudoinverse
+        };
+
+        for (auto algo : algorithms) {
+            SECTION("Solve using algorithm: " + to_string_algorithm(algo)) {
+                Vector weights(3);
+                double rcond = -1.0;
+                std::cout << "------------" << std::endl;
+                std::cout << A << std::endl;
+                std::cout << "expectedWeights: " << expectedWeights << std::endl;
+                ols.solve(A, B, weights, rcond, algo);
+                std::cout << "weights: " << weights << std::endl;
+                Vector B_computed = A * weights;
+                std::cout << "B: " << B<< std::endl;
+                std::cout << "B_computed: " << B_computed << std::endl;
+                vectorsAreClose(weights, expectedWeights, 1e-4);
+            }
+        }
     }
-  };
-
-  // Prepare test data with corrected matrix initialization
-  Matrix A(3, 3, {1, 4, 7,   // First column
-                  2, 5, 8,   // Second column
-                  3, 6, 10}); // Third column
-
-  Vector expected_weights = {1, 2, 3};
-
-  // Recompute B = A * expected_weights
-  Vector B = A * expected_weights;
-
-  SECTION("Solve using default algorithm (GELSD) without workspace") {
-    Vector weights(3);
-    OLS::solve(A, B, weights);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Solve using GELS algorithm without workspace") {
-    Vector weights(3);
-    OLS::solve(A, B, weights, -1.0, OLS::Algorithm::GELS);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Solve using Cholesky algorithm without workspace") {
-    Vector weights(3);
-    OLS::solve(A, B, weights, -1.0, OLS::Algorithm::Cholesky);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Solve using custom SVD implementation without workspace") {
-    Vector weights(3);
-    OLS::solve(A, B, weights, -1.0, OLS::Algorithm::SVD);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Solve using default algorithm (GELSD) with workspace") {
-    Vector weights(3);
-    OLS::Workspace ws;
-    OLS::solve_with_workspace(A, B, weights, ws);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Solve using GELS algorithm with workspace") {
-    Vector weights(3);
-    OLS::Workspace ws;
-    OLS::solve_with_workspace(A, B, weights, ws, -1.0, OLS::Algorithm::GELS);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Solve using Cholesky algorithm with workspace") {
-    Vector weights(3);
-    OLS::Workspace ws;
-    OLS::solve_with_workspace(A, B, weights, ws, -1.0, OLS::Algorithm::Cholesky);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Solve using custom SVD implementation with workspace") {
-    Vector weights(3);
-    OLS::Workspace ws;
-    OLS::solve_with_workspace(A, B, weights, ws, -1.0, OLS::Algorithm::SVD);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Test handling of rank-deficient matrix") {
-    // Corrected initialization in column-major order
-    Matrix A_singular(3, 3, {
-      1, 4, 7,  // First column
-      2, 5, 8,  // Second column
-      3, 6, 9   // Third column
-    });
-
-    // Make second column a multiple of the first: second column = 2 * first column
-    for (int i = 0; i < 3; ++i) {
-      A_singular(i,1) = 2 * A_singular(i,0);
+
+    SECTION("Solve with overdetermined system (m > n)") {
+        int m = 5;
+        int n = 3;
+        Matrix A(m, n);
+
+        // Initialize A (column-major order)
+        A(0,0) = 1;  A(1,0) = 4;  A(2,0) = 7;  A(3,0) = 1;  A(4,0) = 0; // First column
+        A(0,1) = 2;  A(1,1) = 5;  A(2,1) = 8;  A(3,1) = 0;  A(4,1) = 1; // Second column
+        A(0,2) = 3;  A(1,2) = 6;  A(2,2) = 10; A(3,2) = 0;  A(4,2) = 0; // Third column
+
+        Matrix A_cpy = A;
+        Vector expectedWeights = {1, 2, 3};
+
+        // Compute B = A * expectedWeights
+        Vector B = A * expectedWeights;
+        Vector B_cpy = B;
+
+        tadah::models::OLS ols;
+
+        Vector weights(n);
+        ols.solve(A, B, weights);
+
+        // Compute B_computed = A * weights
+        Vector B_computed = A_cpy * weights;
+
+        // Check residual
+        for (int i = 0; i < m; ++i) {
+            REQUIRE(std::abs(B_cpy[i] - B_computed[i]) < 1e-3);
+        }
     }
 
-    // Recompute B_singular to match the modified A_singular and expected weights
-    Vector expected_weights_singular = {1, 0, 3}; // Since second column is dependent, set weight[1] = 0
-    Vector B_singular = A_singular * expected_weights_singular;
+    SECTION("Solve with underdetermined system (m < n)") {
+        int m = 3;
+        int n = 5;
+        Matrix A(m, n);
+
+        // Initialize A (column-major order)
+        A(0,0) = 1;  A(1,0) = 6;   A(2,0) = 11;  // First column
+        A(0,1) = 2;  A(1,1) = 7;   A(2,1) = 12;  // Second column
+        A(0,2) = 3;  A(1,2) = 8;   A(2,2) = 13;  // Third column
+        A(0,3) = 4;  A(1,3) = 9;   A(2,3) = 14;  // Fourth column
+        A(0,4) = 5;  A(1,4) = 10;  A(2,4) = 15;  // Fifth column
 
-    Matrix A_singular2 = A_singular;
-    Vector B_singular2 = B_singular;
+        Matrix A_cpy = A;
+        Vector B = {1, 2, 3};
+        Vector B_cpy = B;
 
-    Vector weights(3);
+        tadah::models::OLS ols;
 
-    // Solve using GELSD, which can handle rank deficiency
-    OLS::solve(A_singular, B_singular, weights);
+        Vector weights(n);
+        ols.solve(A, B, weights);
 
-    // Since the system is rank-deficient, the solution is not unique
-    // We can check that A_singular * weights is close to B_singular
-    Vector B_computed = A_singular2 * weights;
+        // Compute B_computed = A * weights
+        Vector B_computed = A_cpy * weights;
 
-    // Increase tolerance due to potential numerical inaccuracies
-    vectors_are_close(B_computed, B_singular2, 1e-4);
-  }
+        // Check that B_computed matches B
+        vectorsAreClose(B_computed, B_cpy, 1e-4);
+    }
 
-  SECTION("Test solve with overdetermined system") {
-    // Corrected matrix initialization in column-major order
-    int m = 5;
-    int n = 3;
-    Matrix A_over(m, n, {
-      1, 4, 7, 1, 0, // First column
-      2, 5, 8, 0, 1, // Second column
-      3, 6, 10,0, 0  // Third column
-    });
-    Matrix A_over2 = A_over;
+    SECTION("Solve with rank-deficient matrix") {
+        Matrix A(3, 3);
 
-    Vector expected_weights_over = {1, 2, 3};
+        // Initialize A to be rank-deficient (second column is 2 * first column)
+        A(0,0) = 1;  A(1,0) = 2;  A(2,0) = 3; // First column
+        A(0,1) = 2;  A(1,1) = 4;  A(2,1) = 6; // Second column (2 * first column)
+        A(0,2) = 3;  A(1,2) = 6;  A(2,2) = 9; // Third column
 
-    // Recompute B_over
-    Vector B_over = A_over * expected_weights_over;
-    Vector B_over2 = B_over;
+        Matrix A_cpy = A;
+        Vector expectedWeights = {1, 0, 3};
+        Vector B = A * expectedWeights;
+        Vector B_cpy = B;
 
-    Vector weights(n);
+        tadah::models::OLS ols;
 
-    OLS::solve(A_over, B_over, weights);
+        Vector weights(3);
+        ols.solve(A, B, weights);
 
-    // Compute B_computed = A_over * weights
-    Vector B_computed = A_over2 * weights;
+        // Since the system is rank-deficient, the solution is not unique
+        // Verify that A * weights is close to B
 
-    // Check that the residual is small
-    for (int i = 0; i < m; ++i) {
-      REQUIRE(std::abs(B_over2[i] - B_computed[i]) < 1e-3); // Adjusted tolerance
+        Vector B_computed = A_cpy * weights;
+        vectorsAreClose(B_computed, B_cpy, 1e-4);
     }
-  }
-
-  SECTION("Test solve with underdetermined system") {
-    // Corrected matrix initialization in column-major order
-    int m = 3;
-    int n = 5;
-    Matrix A_under(m, n, {
-      1, 6, 11, // First column
-      2, 7, 12, // Second column
-      3, 8, 13, // Third column
-      4, 9, 14, // Fourth column
-      5,10,15  // Fifth column
-    });
-
-    Vector B_under = {1, 2, 3};
-    // Enlarge B_under to size max(m, n) by appending zeros
-    B_under.resize(std::max(m, n));
-
-    Vector weights(n);
-
-    OLS::solve(A_under, B_under, weights);
-
-    // Compute B_computed = A_under * weights
-    Vector B_computed = A_under * weights;
-
-    // Since B_under has been resized, compare only the first m elements
-    Vector B_under_trimmed = Vector(B_under.ptr(), m);
-    Vector B_computed_trimmed = Vector(B_computed.ptr(), m);
-
-    vectors_are_close(B_computed_trimmed, B_under_trimmed, 1e-4);
-  }
-
-  SECTION("Test solve_with_svd without workspace") {
-    // Assuming SVDType is properly defined and we have access to an SVD implementation
-    // Replace 'SVDType' with your actual SVD class or struct
-
-    // For this test, we'll use the custom SVD algorithm on matrix A
-
-    // Placeholder SVDType and methods
-    struct SVDType {
-      double* getU() { return U.data(); }
-      double* getS() { return S.data(); }
-      double* getVT() { return VT.data(); }
-      std::vector<double> U;
-      std::vector<double> S;
-      std::vector<double> VT;
-    };
 
-    // Compute SVD of A
-    SVDType svd;
-    int m = A.rows();
-    int n = A.cols();
-    int min_mn = std::min(m, n);
-    svd.U.resize(m * min_mn);
-    svd.S.resize(min_mn);
-    svd.VT.resize(min_mn * n);
-
-    // Prepare parameters for LAPACK dgesdd_
-    char jobz = 'S';
-    int lda = m;
-    int ldu = m;
-    int ldvt = min_mn;
-    int info;
-    std::vector<double> A_copy(A.ptr(), A.ptr() + m * n);
-
-    // Workspace query
-    int lwork = -1;
-    double wkopt;
-    std::vector<int> iwork(8 * min_mn);
-    ::dgesdd_(&jobz, &m, &n, A_copy.data(), &lda, svd.S.data(), svd.U.data(), &ldu, svd.VT.data(), &ldvt, &wkopt, &lwork, iwork.data(), &info);
-
-    lwork = static_cast<int>(wkopt);
-    std::vector<double> work(lwork);
-
-    // Compute SVD
-    ::dgesdd_(&jobz, &m, &n, A_copy.data(), &lda, svd.S.data(), svd.U.data(), &ldu, svd.VT.data(), &ldvt, work.data(), &lwork, iwork.data(), &info);
-
-    REQUIRE(info == 0);
-
-    Vector weights(3);
-    OLS::solve_with_svd(m, n, B, weights, svd);
-
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Test solve_with_svd with workspace") {
-    // Similar to the previous section but using workspace
-    struct SVDType {
-      double* getU() { return U.data(); }
-      double* getS() { return S.data(); }
-      double* getVT() { return VT.data(); }
-      std::vector<double> U;
-      std::vector<double> S;
-      std::vector<double> VT;
-    };
+    SECTION("Testing solveWithPseudoinverse method") {
+        Matrix A(3, 3);
+        // Initialize A (column-major order)
+        A(0,0) = 1; A(1,0) = 4; A(2,0) = 7;   // First column
+        A(0,1) = 2; A(1,1) = 5; A(2,1) = 8;   // Second column
+        A(0,2) = 3; A(1,2) = 6; A(2,2) = 10;  // Third column
 
-    // Compute SVD of A
-    SVDType svd;
-    int m = A.rows();
-    int n = A.cols();
-    int min_mn = std::min(m, n);
-    svd.U.resize(m * min_mn);
-    svd.S.resize(min_mn);
-    svd.VT.resize(min_mn * n);
-
-    // Prepare parameters for LAPACK dgesdd_
-    char jobz = 'S';
-    int lda = m;
-    int ldu = m;
-    int ldvt = min_mn;
-    int info;
-    std::vector<double> A_copy(A.ptr(), A.ptr() + m * n);
-
-    // Workspace query
-    int lwork = -1;
-    double wkopt;
-    std::vector<int> iwork(8 * min_mn);
-    ::dgesdd_(&jobz, &m, &n, A_copy.data(), &lda, svd.S.data(), svd.U.data(), &ldu, svd.VT.data(), &ldvt, &wkopt, &lwork, iwork.data(), &info);
-
-    lwork = static_cast<int>(wkopt);
-    std::vector<double> work(lwork);
-
-    // Compute SVD
-    ::dgesdd_(&jobz, &m, &n, A_copy.data(), &lda, svd.S.data(), svd.U.data(), &ldu, svd.VT.data(), &ldvt, work.data(), &lwork, iwork.data(), &info);
-
-    REQUIRE(info == 0);
-
-    Vector weights(3);
-    OLS::Workspace ws;
-    double rcond = -1.0;
-    double lambda = 0.0;
-
-    OLS::solve_with_svd_with_workspace(m, n, B, weights, svd, ws, rcond, lambda);
-
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
-
-  SECTION("Test solve_with_svd with regularization (lambda)") {
-    // Using a non-zero lambda for regularization
-    struct SVDType {
-      double* getU() { return U.data(); }
-      double* getS() { return S.data(); }
-      double* getVT() { return VT.data(); }
-      std::vector<double> U;
-      std::vector<double> S;
-      std::vector<double> VT;
-    };
+        Vector expectedWeights = {1, 2, 3};
+
+        // Compute B = A * expectedWeights
+        Vector B = A * expectedWeights;
+
+        tadah::models::OLS ols;
+        tadah::models::SVD svd;
 
-    // Compute SVD of A
-    SVDType svd;
-    int m = A.rows();
-    int n = A.cols();
-    int min_mn = std::min(m, n);
-    svd.U.resize(m * min_mn);
-    svd.S.resize(min_mn);
-    svd.VT.resize(min_mn * n);
+        // Compute SVD of A
+        svd.compute(A);
 
-    // Prepare parameters for LAPACK dgesdd_
-    char jobz = 'S';
-    int lda = m;
-    int ldu = m;
-    int ldvt = min_mn;
-    int info;
-    std::vector<double> A_copy(A.ptr(), A.ptr() + m * n);
+        Vector weights(3);
+        ols.solveWithPseudoinverse(A.rows(), A.cols(), B, weights, svd.getSVDWorkspace());
 
-    // Workspace query
-    int lwork = -1;
-    double wkopt;
-    std::vector<int> iwork(8 * min_mn);
-    ::dgesdd_(&jobz, &m, &n, A_copy.data(), &lda, svd.S.data(), svd.U.data(), &ldu, svd.VT.data(), &ldvt, &wkopt, &lwork, iwork.data(), &info);
+        vectorsAreClose(weights, expectedWeights, 1e-4);
+    }
+
+    SECTION("Testing solveWithPseudoinverse with regularization") {
+        Matrix A(3, 3);
+        // Initialize A (column-major order)
+        A(0,0) = 1; A(1,0) = 4; A(2,0) = 7;   // First column
+        A(0,1) = 2; A(1,1) = 5; A(2,1) = 8;   // Second column
+        A(0,2) = 3; A(1,2) = 6; A(2,2) = 10;  // Third column
 
-    lwork = static_cast<int>(wkopt);
-    std::vector<double> work(lwork);
+        Matrix A_cpy = A;
+        Vector expectedWeights = {1, 2, 3};
 
-    // Compute SVD
-    ::dgesdd_(&jobz, &m, &n, A_copy.data(), &lda, svd.S.data(), svd.U.data(), &ldu, svd.VT.data(), &ldvt, work.data(), &lwork, iwork.data(), &info);
+        // Compute B = A * expectedWeights
+        Vector B = A * expectedWeights;
+        Vector B_cpy = B;
 
-    REQUIRE(info == 0);
+        tadah::models::OLS ols;
+        tadah::models::SVD svd;
 
-    Vector weights(3);
-    double rcond = -1.0;
-    double lambda = 0.0005; // Regularization parameter
+        // Compute SVD of A
+        svd.compute(A);
 
-    OLS::solve_with_svd(m, n, B, weights, svd, rcond, lambda);
+        Vector weights(3);
+        double rcond = -1.0;
+        double lambda = 0.0001; // Regularization parameter
 
-    // Since lambda introduces regularization, weights may not match exactly
-    // but for this test, we can still check that the solution is reasonable
+        ols.solveWithPseudoinverse(A.rows(), A.cols(), B, weights, svd.getSVDWorkspace(), rcond, lambda);
 
-    Vector B_computed = A * weights;
+        // Since lambda introduces regularization, the weights may not match exactly
+        // Compute residual
+        Vector B_computed = A_cpy * weights;
+        double residualNorm = 0.0;
 
-    // Compute residual norm
-    double residual_norm = 0.0;
-    for (int i = 0; i < m; ++i) {
-      residual_norm += std::pow(B[i] - B_computed[i], 2);
+        for (size_t i = 0; i < A.rows(); ++i) {
+            residualNorm += std::pow(B_cpy[i] - B_computed[i], 2);
+        }
+        residualNorm = std::sqrt(residualNorm);
+
+        REQUIRE(residualNorm < 1e-3);
     }
-    residual_norm = std::sqrt(residual_norm);
 
-    REQUIRE(residual_norm < 1e-3); // Adjust tolerance as needed
-  }
+    SECTION("Testing with workspace manager") {
+        Matrix A(3, 3);
+        // Initialize A (column-major order)
+        A(0,0) = 1; A(1,0) = 4; A(2,0) = 7;   // First column
+        A(0,1) = 2; A(1,1) = 5; A(2,1) = 8;   // Second column
+        A(0,2) = 3; A(1,2) = 6; A(2,2) = 10;  // Third column
+
+        Vector expectedWeights = {1, 2, 3};
+        Vector B = A * expectedWeights;
 
-  SECTION("Test solve with Cholesky algorithm on ill-conditioned matrix") {
-    // Create an ill-conditioned matrix
-    Matrix A_ill(3, 3, {1e-11, 0,   0,
-                        0,   1e-11, 0,
-                        0,    0,   1e-11});
+        // Create a workspace manager
+        tadah::models::memory::ModelsWorkspaceManager workspaceManager;
 
-    Vector B_ill = {1e-11, 1e-11, 1e-11};
-    Vector weights(3);
+        // Create OLS instance with the workspace manager
+        tadah::models::OLS ols(workspaceManager);
+
+        Vector weights(3);
+        ols.solve(A, B, weights);
+
+        vectorsAreClose(weights, expectedWeights, 1e-4);
+
+        // Reuse the workspace manager for another computation
+        Matrix A2(3, 3);
+        A2(0,0) = 2; A2(1,0) = 5; A2(2,0) = 8;
+        A2(0,1) = 3; A2(1,1) = 6; A2(2,1) = 9;
+        A2(0,2) = 4; A2(1,2) = 7; A2(2,2) = 11;
+
+        Vector expectedWeights2 = {2, 3, 4};
+        Vector B2 = A2 * expectedWeights2;
+
+        Vector weights2(3);
+        ols.solve(A2, B2, weights2);
+
+        vectorsAreClose(weights2, expectedWeights2, 1e-4);
+    }
 
-    // Solve using Cholesky algorithm
-    OLS::solve(A_ill, B_ill, weights, -1.0, OLS::Algorithm::Cholesky);
-    Vector B_computed = A_ill * weights;
-    vectors_are_close(B_ill, B_computed, 1e-11);
-  }
+    SECTION("Testing exception handling with invalid inputs (NormalEquations)") {
+        Matrix A(3, 3);
+        // Initialize A (column-major order) with zeros to create a singular matrix
+        for (int i = 0; i < 3; ++i) {
+            for (int j = 0; j < 3; ++j) {
+                A(i,j) = 0.0;
+            }
+        }
 
-  SECTION("Test solve with Cholesky algorithm on well-conditioned matrix") {
-    Matrix A_well = A; // Use the original well-conditioned matrix
-    Vector B_well = B;
-    Vector weights(3);
+        Vector B = {0, 0, 0};
+        Vector weights(3);
 
-    OLS::solve(A_well, B_well, weights, -1.0, OLS::Algorithm::Cholesky);
-    vectors_are_close(weights, expected_weights, 1e-4);
-  }
+        tadah::models::OLS ols;
 
-  // Rest of the tests...
+        // Expect an exception due to singular A^T A matrix in NormalEquations algorithm
+        REQUIRE_THROWS_AS(ols.solve(A, B, weights, -1.0, tadah::models::Algorithm::NormalEquations), std::runtime_error);
+    }
+
+    SECTION("Testing with ill-conditioned matrix") {
+        Matrix A(3, 3);
+
+        // Create an ill-conditioned matrix
+        A(0,0) = 1e-10; A(1,0) = 0;      A(2,0) = 0;
+        A(0,1) = 0;      A(1,1) = 1e-10; A(2,1) = 0;
+        A(0,2) = 0;      A(1,2) = 0;      A(2,2) = 1e-10;
+
+        Vector expectedWeights = {1, 2, 3};
+        Vector B = A * expectedWeights;
+
+        tadah::models::OLS ols;
+
+        Vector weights(3);
+        ols.solve(A, B, weights, -1.0, tadah::models::Algorithm::GELSD);
+
+        vectorsAreClose(weights, expectedWeights, 1e-2);
+    }
+
+    SECTION("Testing NormalEquations algorithm with rank-deficient matrix") {
+        Matrix A(2, 2);
+
+        // Initialize A (column-major order) with a rank-deficient matrix
+        A(0,0) = 1;  A(1,0) = 2;  // First column
+        A(0,1) = 2;  A(1,1) = 4;  // Second column (2 * first column)
+
+        Vector B = {5, 10}; // Compute B = A * [1; 2], but it's arbitrary here
+        Vector weights(2);
+
+        tadah::models::OLS ols;
+
+        // Expect an exception due to singular A^T A matrix
+        REQUIRE_THROWS_AS(ols.solve(A, B, weights, -1.0, tadah::models::Algorithm::NormalEquations), std::runtime_error);
+    }
+
+    SECTION("Testing solve with minimum input") {
+        Matrix A(1, 1);
+        A(0,0) = 2;
+
+        Vector B = {4};
+        Vector weights(1);
+
+        tadah::models::OLS ols;
+
+        ols.solve(A, B, weights);
+
+        REQUIRE(weights[0] == Approx(2.0).epsilon(1e-6));
+    }
+
+    SECTION("Testing solve with large matrix") {
+        int m = 100;
+        int n = 50;
+        Matrix A(m, n);
+
+        // Initialize A with random values
+        for (int j = 0; j < n; ++j) {
+            for (int i = 0; i < m; ++i) {
+                A(i,j) = static_cast<double>(std::rand()) / RAND_MAX;
+            }
+        }
+
+        Matrix A_cpy = A;
+
+        Vector expectedWeights(n);
+        for (int i = 0; i < n; ++i) {
+            expectedWeights[i] = static_cast<double>(std::rand()) / RAND_MAX;
+        }
+
+        Vector B = A * expectedWeights;
+        Vector B_cpy = B;
+
+        tadah::models::OLS ols;
+
+        Vector weights(n);
+        ols.solve(A, B, weights);
+
+        // Compute residual norm
+        Vector B_computed = A_cpy * weights;
+        double residualNorm = 0.0;
+        for (int i = 0; i < m; ++i) {
+            residualNorm += std::pow(B_cpy[i] - B_computed[i], 2);
+        }
+        residualNorm = std::sqrt(residualNorm);
+
+        REQUIRE(residualNorm < 1e-2);
+    }
 }
 
diff --git a/tests/test_ridge_regression.cpp b/tests/test_ridge_regression.cpp
deleted file mode 100644
index 064a0cd..0000000
--- a/tests/test_ridge_regression.cpp
+++ /dev/null
@@ -1,31 +0,0 @@
-#include "catch2/catch.hpp"
-#include <tadah/core/maths.h>
-#include <tadah/models/ridge_regression.h>
-#include <iomanip>
-
-TEST_CASE("Testing RR") {
-  Matrix Phi1(3, 3, {1, 4, 7, 2, 5, 8, 3, 6, 9});
-  Matrix Phi2(3, 3, {1, 4, 7, 2, 5, 8, 3, 6, 9});
-  aed_type b1(3);
-  b1[0]=1;
-  b1[1]=1;
-  b1[2]=1;
-  aed_type w(3);
-  aed_type b2=b1;
-
-  SECTION("lambda zero") {
-    double lambda = 0;
-    RidgeRegression::solve(Phi1,b1,w, lambda, 1e-8);
-    aed_type p= Phi2*w;
-    REQUIRE_THAT(p[0], Catch::Matchers::WithinRel(b2[0]));
-    REQUIRE_THAT(p[1], Catch::Matchers::WithinRel(b2[1]));
-    REQUIRE_THAT(p[2], Catch::Matchers::WithinRel(b2[2]));
-  }
-  SECTION("small lambda ") {
-    double lambda = 1e-10;
-    RidgeRegression::solve(Phi1,b1,w, lambda, 1e-8);
-    aed_type p= Phi2*w;
-    REQUIRE(p.isApprox(b2,lambda));
-  }
-}
-
diff --git a/tests/test_svd.cpp b/tests/test_svd.cpp
index d434096..e968aee 100644
--- a/tests/test_svd.cpp
+++ b/tests/test_svd.cpp
@@ -1,98 +1,228 @@
+// test_svd.cpp
+
 #include "catch2/catch.hpp"
+#include <algorithm>
+#include <cmath>
 #include <tadah/core/maths.h>
+#include <tadah/models/memory/ModelsWorkspaceManager.h>
 #include <tadah/models/svd.h>
 
-void multiplyMatrices(const double* A, const double* B, double* C, int m, int n, int p) {
-  // C = A * B
-  // A is m x n, B is n x p, C is m x p
-  for (int i = 0; i < m; ++i) {
-    for (int j = 0; j < p; ++j) {
-      C[i + j * m] = 0;
-      for (int k = 0; k < n; ++k) {
-        C[i + j * m] += A[i + k * m] * B[k + j * n]; // Column-major order
-      }
+TEST_CASE("Testing SVD") {
+  using Matrix = Matrix_T<>;
+  /*using Vector = aed_type;*/
+
+  // Test data
+
+  SECTION("SVD of a rectangular matrix (m > n)") {
+    int m = 5;
+    int n = 3;
+    Matrix A(m, n);
+    // Initialize A with arbitrary values
+    A(0, 0) = 1;
+    A(0, 1) = 2;
+    A(0, 2) = 3;
+    A(1, 0) = 4;
+    A(1, 1) = 5;
+    A(1, 2) = 6;
+    A(2, 0) = 7;
+    A(2, 1) = 8;
+    A(2, 2) = 9;
+    A(3, 0) = 10;
+    A(3, 1) = 11;
+    A(3, 2) = 12;
+    A(4, 0) = 13;
+    A(4, 1) = 14;
+    A(4, 2) = 15;
+
+    int minMN = std::min(m, n);
+
+    tadah::models::SVD svd;
+    svd.compute(A);
+
+    // For this test, we can check that the singular values are non-negative and
+    // in descending order
+    double *S_computed = svd.getS();
+    for (int i = 0; i < minMN - 1; ++i) {
+      REQUIRE(S_computed[i] >= S_computed[i + 1]);
+      REQUIRE(S_computed[i] >= 0);
     }
   }
-}
-
 
-TEST_CASE("SVD Constructor and Destruction") {
-  Matrix Phi(3, 3, {1, 2, 3, 4, 5, 6, 7, 8, 9});
-  SVD svd(Phi);
+  SECTION("SVD of a rectangular matrix (m < n)") {
+    int m = 3;
+    int n = 5;
+    Matrix A(m, n);
+    // Initialize A with arbitrary values
+    A(0, 0) = 1;
+    A(0, 1) = 2;
+    A(0, 2) = 3;
+    A(0, 3) = 4;
+    A(0, 4) = 5;
+    A(1, 0) = 6;
+    A(1, 1) = 7;
+    A(1, 2) = 8;
+    A(1, 3) = 9;
+    A(1, 4) = 10;
+    A(2, 0) = 11;
+    A(2, 1) = 12;
+    A(2, 2) = 13;
+    A(2, 3) = 14;
+    A(2, 4) = 15;
+
+    int minMN = std::min(m, n);
+
+    tadah::models::SVD svd;
+    svd.compute(A);
+
+    // Check that singular values are non-negative and in descending order
+    double *S_computed = svd.getS();
+    for (int i = 0; i < minMN - 1; ++i) {
+      REQUIRE(S_computed[i] >= S_computed[i + 1]);
+      REQUIRE(S_computed[i] >= 0);
+    }
+  }
 
-  REQUIRE(svd.get_info() == 0);
-  REQUIRE(svd.sizeU() == 9);
-  REQUIRE(svd.sizeVT() == 9);
+  SECTION("SVD of a rank-deficient matrix") {
+    Matrix A(4, 4);
+    // Create a rank-deficient matrix (rank 2)
+    A(0, 0) = 1;
+    A(0, 1) = 2;
+    A(0, 2) = 3;
+    A(0, 3) = 4;
+    A(1, 0) = 2;
+    A(1, 1) = 4;
+    A(1, 2) = 6;
+    A(1, 3) = 8;
+    A(2, 0) = 3;
+    A(2, 1) = 6;
+    A(2, 2) = 9;
+    A(2, 3) = 12;
+    A(3, 0) = 4;
+    A(3, 1) = 8;
+    A(3, 2) = 12;
+    A(3, 3) = 16;
+
+    tadah::models::SVD svd;
+    svd.compute(A);
+
+    // Singular values should have two zeros (or very close to zero)
+    double *S_computed = svd.getS();
+    int minMN = std::min(A.rows(), A.cols());
+
+    int zero_count = 0;
+    double tolerance = 1e-6;
+    for (int i = 0; i < minMN; ++i) {
+      if (S_computed[i] < tolerance) {
+        zero_count++;
+      }
+    }
 
-  // Ensure that singular values are correctly allocated
-  double* s = svd.getS();
-  for (size_t i = 0; i < svd.sizeS(); ++i) {
-    REQUIRE(s[i] >= 0.0);
+    REQUIRE(zero_count >=
+            2); // At least two singular values should be near zero
   }
-}
 
-TEST_CASE("SVD Copy Constructor") {
-  Matrix Phi(3, 3, {1, 2, 3, 4, 5, 6, 7, 8, 9});
-  SVD svd(Phi);    
-  SVD svd_copy(svd);
-
-  REQUIRE(svd_copy.get_info() == svd.get_info());
-  REQUIRE(svd_copy.sizeU() == svd.sizeU());
-  REQUIRE(svd_copy.sizeVT() == svd.sizeVT());
+  SECTION("SVD of an identity matrix") {
+    int size = 5;
+    Matrix A(size, size);
+    // Initialize A as an identity matrix
+    for (int i = 0; i < size; ++i) {
+      for (int j = 0; j < size; ++j) {
+        A(i, j) = (i == j) ? 1.0 : 0.0;
+      }
+    }
 
-  // Verify deep copy
-  for (size_t i = 0; i < svd.sizeS(); ++i) {
-    REQUIRE(svd_copy.getS()[i] == svd.getS()[i]);
-  }
-}
+    tadah::models::SVD svd;
+    svd.compute(A);
 
-TEST_CASE("SVD Numeric Accuracy 2x2") {
-  Matrix Phi(2, 2, {1, 2, 3, 4});
-  SVD svd(Phi);
-  REQUIRE(svd.get_info() == 0);
+    double *S_computed = svd.getS();
 
-  double* S = svd.getS();
+    // All singular values should be 1
+    for (int i = 0; i < size; ++i) {
+      REQUIRE(std::abs(S_computed[i] - 1.0) < 1e-6);
+    }
 
-  // S are sorted S[0] > S[1]
-  REQUIRE_THAT(5.464985704219043, Catch::Matchers::WithinRel(S[0]));
-  REQUIRE_THAT(0.365966190626258, Catch::Matchers::WithinRel(S[1]));
-}
-TEST_CASE("SVD Numeric Accuracy - Column-Major Order") {
-  Matrix Phi(3, 3, {1, 2, 3, 4, 5, 6, 7, 8, 9});
-  Matrix Phi_cpy(Phi); // as svd destroys original matrix
-  SVD svd(Phi);
+    // U and VT should be identity matrices (up to sign)
+    double *U = svd.getU();
+    double *VT = svd.getVT();
 
+    // Verify that U is approximately an identity matrix
+    for (int i = 0; i < size; ++i) {
+      for (int j = 0; j < size; ++j) {
+        double expected = (i == j) ? 1.0 : 0.0;
+        REQUIRE(std::abs(std::abs(U[i * size + j]) - expected) < 1e-6);
+      }
+    }
 
-  REQUIRE(svd.get_info() == 0);
+    // Verify that VT is approximately an identity matrix
+    for (int i = 0; i < size; ++i) {
+      for (int j = 0; j < size; ++j) {
+        double expected = (i == j) ? 1.0 : 0.0;
+        REQUIRE(std::abs(std::abs(VT[i * size + j]) - expected) < 1e-6);
+      }
+    }
+  }
 
-  double* U = svd.getU();
-  double* S = svd.getS();
-  double* VT = svd.getVT();
+  SECTION("Access to SVDWorkspace") {
+    Matrix A(3, 3);
+    A(0, 0) = 1;
+    A(0, 1) = 0;
+    A(0, 2) = 0;
+    A(1, 0) = 0;
+    A(1, 1) = 2;
+    A(1, 2) = 0;
+    A(2, 0) = 0;
+    A(2, 1) = 0;
+    A(2, 2) = 3;
+
+    tadah::models::SVD svd;
+    svd.compute(A);
+
+    tadah::models::memory::SVDWorkspace &svdWorkspace = svd.getSVDWorkspace();
+
+    // Verify that the workspace dimensions match
+    REQUIRE(svdWorkspace.m == 3);
+    REQUIRE(svdWorkspace.n == 3);
+
+    // Verify that the singular values match the expected values
+    double *S_computed = svd.getS();
+    REQUIRE(std::abs(S_computed[0] - 3.0) < 1e-6);
+    REQUIRE(std::abs(S_computed[1] - 2.0) < 1e-6);
+    REQUIRE(std::abs(S_computed[2] - 1.0) < 1e-6);
+  }
 
-  int m = svd.shapeU().first;
-  int n = svd.shapeU().second;
+  SECTION("Test with workspace manager") {
+    // Create a workspace manager
+    tadah::models::memory::ModelsWorkspaceManager workspaceManager;
 
-  Matrix mU (U, m,n);
-  Matrix mVT (VT,n,n);
+    // Create an SVD instance with the workspace manager
+    tadah::models::SVD svd(workspaceManager);
 
-  aed_type vec(S, n);
-  REQUIRE_THAT(16.848103352614209, Catch::Matchers::WithinRel(S[0]));
-  REQUIRE_THAT(1.068369514554709, Catch::Matchers::WithinRel(S[1]));
-  REQUIRE_THAT(0, Catch::Matchers::WithinAbs(S[2],1e-15));
+    Matrix A(4, 4);
+    // Initialize A with arbitrary values
+    for (size_t i = 0; i < A.rows(); ++i) {
+      for (size_t j = 0; j < A.cols(); ++j) {
+        A(i, j) = static_cast<double>(i + j + 1);
+      }
+    }
 
-  double USV[9];
-  double SV[9];
+    svd.compute(A);
 
-  // S * VT: We can multiply S into VT directly since S is diagonal
-  for (int i = 0; i < 3; ++i) {
-    for (int j = 0; j < 3; ++j) {
-      SV[i + j * 3] = S[i] * VT[i + j * 3]; // Diagonal multiplication
+    // Perform another computation to test workspace reuse
+    Matrix B(4, 4);
+    for (size_t i = 0; i < B.rows(); ++i) {
+      for (size_t j = 0; j < B.cols(); ++j) {
+        B(i, j) = static_cast<double>((i + 1) * (j + 1));
+      }
     }
-  }
 
-  multiplyMatrices(U, SV, USV, 3, 3, 3);
+    svd.compute(B);
 
-  for (int i = 0; i < 9; ++i)
-    REQUIRE_THAT(USV[i], Catch::Matchers::WithinRel(Phi_cpy[i]));
+    // If no errors occur, the test passes
+    REQUIRE(true);
+  }
 
+  // Additional tests can be added to cover more cases, such as:
+  // - Matrices with complex entries (if supported)
+  // - Very large matrices to test performance and memory management
+  // - Matrices with known singular values to validate accuracy
 }
-- 
GitLab