diff --git a/include/tadah/models/linear_regressor.h b/include/tadah/models/linear_regressor.h
index 0485982bb5bb8965ad0167b25c420a3124ae9fcf..92c47ea35cc06a7c6b46e3e1caed9f1a83887e28 100644
--- a/include/tadah/models/linear_regressor.h
+++ b/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;
diff --git a/include/tadah/models/ols.h b/include/tadah/models/ols.h
index fce36561237f9514da52294b144d77fabf2ef2d5..e6c11f6650606439928f52d6c9816e4baa179aaa 100644
--- a/include/tadah/models/ols.h
+++ b/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) {
diff --git a/src/ols.cpp b/src/ols.cpp
index 22836f1582fc49dfac1cf6f6a08669dd28d0f0a4..cd2eced12464f122cb6f6442eb415b5eef3cbbaf 100644
--- a/src/ols.cpp
+++ b/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;
}
-