Not working SVD implementation in OLS, compiles but produces NaN

include/tadah/models/linear_regressor.h

@@ -65,7 +65,7 @@ public:
     }
     else if (oalgo==3) {
       SVD svd = SVD(Phi);;
-      //OLS::Workspace ws;
+      // OLS::Workspace ws;
       if (lambda < 0) {
         double alpha = config.get<double>("ALPHA");
         double beta = config.get<double>("BETA");
@@ -78,16 +78,19 @@ public:
         config.add("BETA", beta);
       }
 
-      //ws.U = svd.getU();
-      //ws.S = svd.getS();
-      //ws.VT = svd.getVT();
 
       Sigma = get_sigma(svd, lambda);
       //std::cout << "Sigma: " << Sigma << std::endl;
       config_add_sigma(config, Sigma);
 
       if (verbose) std::cout << std::endl << "REG LAMBDA: " << lambda << std::endl;
+      // ws.U = svd.getU();
+      // ws.S = svd.getS();
+      // ws.VT = svd.getVT();
+
+      // std::cout << "own_U: " << ws.owns_U << std::endl;
       //OLS::solve_with_workspace(Phi, T, weights, ws, rcond, OLS::Algorithm::SVD);
+      // OLS::solve(Phi, T, weights, rcond, OLS::Algorithm::SVD);
       RidgeRegression::solve(svd, T, weights, lambda, rcond);
       //Sigma = get_sigma(svd, lambda);
       //std::cout << "Sigma: " << Sigma << std::endl;

include/tadah/models/ols.h

@@ -73,6 +73,12 @@ public:
     double* U;    // U matrix from SVD (m x n)
     double* S;    // Singular values array (size min(m,n))
     double* VT;   // V^T matrix from SVD (n x n)
+    int* iwork_svd; // Integer workspace array for SVD
+
+    // Ownership flags for U, S, VT
+    bool owns_U;
+    bool owns_S;
+    bool owns_VT;
 
     // For Cholesky algorithm
     double* AtA;  // A^T * A matrix (n x n)
@@ -179,14 +185,75 @@ void OLS::solve_with_workspace(M& A, V& B, V& weights, Workspace& ws, double rco
       weights[i] = b[i];
     }
   } else if (algorithm == Algorithm::SVD) {
+
     // Custom SVD-based solution
-    // Implementation can be provided or left as an exercise
-    //throw std::runtime_error("SVD algorithm is not implemented in this version.");
-    // Allocate memory for intermediate calculations
+    int min_mn = std::min(m, n);
+
+    // Check if U, S, VT are provided
+    bool have_svd_matrices = (ws.U != nullptr) && (ws.S != nullptr) && (ws.VT != nullptr);
+
+    if (!have_svd_matrices) {
+      // 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;
 
+      // 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));
+      }
+    }
+
+    // Now compute x = V * S⁻¹ * Uᵗ * b
+
+    // Compute Uᵗ * b
+    std::vector<double> U_t_b(min_mn, 0.0);
+    char trans = 'T'; // Transpose U
+    double alpha = 1.0;
+    double beta = 0.0;
+    int incx = 1;
+    int incy = 1;
+    int m_u = m;
+    int n_u = min_mn;
+    ::dgemv_(&trans, &m_u, &n_u, &alpha, ws.U, &m_u, B.ptr(), &incx, &beta, U_t_b.data(), &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;
+
+    for (int i = 0; i < min_mn; ++i) {
+      if (ws.S[i] > thresh) {
+        U_t_b[i] /= ws.S[i];
+        //U_t_b[i] /= ws.S[i] / (ws.S[i] * ws.S[i] + ws.lambda);
+      } else {
+        // Singular value is too small; set result to zero
+        U_t_b[i] = 0.0;
+      }
+    }
+
 
+    // Compute x = V * (Uᵗ * b / S)
+    trans = 'T'; // Transpose VT to get V
+    int m_vt = min_mn;
+    int n_vt = n;
+    ::dgemv_(&trans, &m_vt, &n_vt, &alpha, ws.VT, &m_vt, U_t_b.data(), &incx, &beta, weights.ptr(), &incy);
 
 
   } else if (algorithm == Algorithm::Cholesky) {

src/ols.cpp

@@ -9,8 +9,10 @@
 // Constructor
 OLS::Workspace::Workspace()
   : work(nullptr), lwork(-1), m(0), n(0), algo(Algorithm::GELSD),
-  s(nullptr), iwork(nullptr), U(nullptr), S(nullptr), VT(nullptr),
-  AtA(nullptr), Atb(nullptr) {
+    s(nullptr), iwork(nullptr),
+    U(nullptr), S(nullptr), VT(nullptr), iwork_svd(nullptr),
+    owns_U(false), owns_S(false), owns_VT(false),
+    AtA(nullptr), Atb(nullptr) {
 }
 
 // Destructor
@@ -63,48 +65,38 @@ void OLS::Workspace::allocate(int m_, int n_, Algorithm algorithm) {
     lwork = static_cast<int>(work_query);
     work = new double[lwork];
   } else if (algorithm == Algorithm::SVD) {
-    // Allocate workspace for custom SVD implementation using DGESDD
-    int minmn = std::min(m, n);
-    int maxmn = std::max(m, n);
-
-    // Set jobz parameter to compute singular vectors
-    char jobz = 'S'; // Compute the first min(m,n) singular vectors
-
-    // Leading dimensions
-    int lda = m;      // Leading dimension of A
-    int ldu = m;      // Leading dimension of U
-    int ldvt = minmn; // Leading dimension of VT
-
-    // Allocate arrays for U, S, and VT
-    U = new double[ldu * minmn];       // U is m x min(m,n)
-    S = new double[minmn];             // Singular values vector
-    VT = new double[ldvt * n];         // VT is min(m,n) x n
-
-    // Allocate integer workspace for DGESDD
-    int* iwork = new int[8 * minmn];
-
-    // Workspace query for DGESDD
-    double work_query;
-    int lwork = -1;
-    int info;
 
-    // Create a dummy A matrix for the workspace query
+    // Allocate U, S, VT, and iwork_svd
+    int min_mn = std::min(m, n);
+
+    U = new double[m * min_mn];      // U is m x min(m,n)
+    S = new double[min_mn];          // Singular values
+    VT = new double[min_mn * n];     // VT is min(m,n) x n
+    iwork_svd = new int[8 * min_mn]; // Workspace for dgesdd_
+
+    owns_U = true;
+    owns_S = true;
+    owns_VT = true;
+
+    // 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);
 
-    // Perform the workspace query
-    ::dgesdd_(&jobz, &m, &n, a_dummy, &lda, S, U, &ldu, VT, &ldvt,
-              &work_query, &lwork, iwork, &info);
-
-    // Check for errors
-    if (info != 0) {
-      delete[] iwork;
-      throw std::runtime_error("Error in DGESDD workspace query: INFO = " + std::to_string(info));
+    if (info_query != 0) {
+      throw std::runtime_error("Error in DGESDD (workspace query): info = " + std::to_string(info_query));
     }
 
-    // Allocate the optimal workspace
-    lwork = static_cast<int>(work_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];
@@ -120,21 +112,32 @@ void OLS::Workspace::allocate(int m_, int n_, Algorithm algorithm) {
 
 // release method
 void OLS::Workspace::release() {
+  if (owns_U) {
+    delete[] U;
+    U = nullptr;
+    owns_U = false;
+  }
+  if (owns_S) {
+    delete[] S;
+    S = nullptr;
+    owns_S = false;
+  }
+  if (owns_VT) {
+    delete[] VT;
+    VT = nullptr;
+    owns_VT = false;
+  }
   delete[] s;
   delete[] iwork;
+  delete[] iwork_svd;
   delete[] work;
-  delete[] U;
-  delete[] S;
-  delete[] VT;
   delete[] AtA;
   delete[] Atb;
 
   s = nullptr;
   iwork = nullptr;
+  iwork_svd = nullptr;
   work = nullptr;
-  U = nullptr;
-  S = nullptr;
-  VT = nullptr;
   AtA = nullptr;
   Atb = nullptr;
 
@@ -152,4 +155,3 @@ bool OLS::Workspace::is_sufficient(int m_, int n_, Algorithm algorithm) const {
   // Additional checks can be added here if necessary
   return true;
 }
-