diff --git a/include/tadah/models/linear_regressor.h b/include/tadah/models/linear_regressor.h
index 4e7355fb66fcd165102a4f57cafb63abf2eb432d..54df74aa020a78b7bef9ffff5a60070d85d22061 100644
--- a/include/tadah/models/linear_regressor.h
+++ b/include/tadah/models/linear_regressor.h
@@ -45,22 +45,18 @@ public:
throw std::runtime_error("LAMBDA -1 requires OALGO 3");
}
- if (oalgo==1) {
+ if (lambda > 0 && (oalgo==1 || oalgo==2)) {
+ // Augment the design matrix with the regularization term
size_t j=0;
- if (lambda>0) {
- for (size_t i=Phi.rows()-Phi.cols(); i<Phi.rows();++i) {
- Phi(i,j++) = std::sqrt(lambda);
- }
+ 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) {
- size_t j=0;
- if (lambda>0) {
- for (size_t i=Phi.rows()-Phi.cols(); i<Phi.rows();++i) {
- Phi(i,j++) = std::sqrt(lambda);
- }
- }
OLS::solve(Phi, T, weights, rcond, OLS::Algorithm::GELS);
}
else if (oalgo==3) {
@@ -83,6 +79,9 @@ public:
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 {
throw std::runtime_error("Unsupported OALGO: " + std::to_string(oalgo));
}
diff --git a/include/tadah/models/ols.h b/include/tadah/models/ols.h
index 3c4ec409f24cb6d13ff5d9e8677291dba3ffbda7..1c427a97227a46aefa9e6ac597972b09268f3717 100644
--- a/include/tadah/models/ols.h
+++ b/include/tadah/models/ols.h
@@ -20,7 +20,7 @@
* - **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 using Cholesky decomposition, efficient but may be less stable for ill-conditioned systems.
+ * - **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.
@@ -253,7 +253,8 @@ void OLS::solve_with_workspace(M& A, V& B, V& weights, Workspace& ws, double rco
// Custom SVD-based solution
int min_mn = std::min(m, n);
- if (!ws.U || !ws.S || !ws.VT) {
+ //if (!ws.U || !ws.S || !ws.VT) {
+ if (ws.owns_U || ws.owns_S || ws.owns_VT) {
// Perform SVD
// Ensure workspace is allocated
@@ -273,9 +274,6 @@ void OLS::solve_with_workspace(M& A, V& B, V& weights, Workspace& ws, double rco
throw std::runtime_error("Error in DGESDD: info = " + std::to_string(info));
}
- ws.owns_U = true;
- ws.owns_S = true;
- ws.owns_VT = true;
}
if (ws.U_t_b == nullptr) {
@@ -318,109 +316,141 @@ void OLS::solve_with_workspace(M& A, V& B, V& weights, Workspace& ws, double rco
} else if (algorithm == Algorithm::Cholesky) {
// Cholesky-based solution
- throw std::runtime_error("Cholesky algorithm is not implemented in this version.");
+
+ // 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.");
}
}
-// New methods: solve_with_svd
-
// 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);
+ 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);
+ 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);
+ 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.");
+ 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];
}
+ }
- // 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.");
- }
+ // Threshold for singular values
+ double thresh = rcond > 0.0 ? rcond * smax : 0.0;
- // 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;
-
- // 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;
- }
+ // 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;
}
+ }
- // 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);
+ // 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 // OLS_H
diff --git a/src/ols.cpp b/src/ols.cpp
index 6de4092eaa721bd169b24b9f49f75bc3e9e82956..52bc49f09da67fdf4339aff99d56fa93e4c0b326 100644
--- a/src/ols.cpp
+++ b/src/ols.cpp
@@ -8,156 +8,158 @@
// 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) {
+ : 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();
+ 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);
-
- // Only allocate if we own the data
- if (owns_U) {
- U = new double[m * min_mn]; // U is m x min(m,n)
- }
- if (owns_S) {
- S = new double[min_mn]; // Singular values
- }
- if (owns_VT) {
- VT = new double[min_mn * n]; // VT is min(m,n) x n
- }
- 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(); // 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;
}
-}
-
-// release method
-void OLS::Workspace::release() {
- if (owns_U && U != nullptr) {
- delete[] U;
- U = nullptr;
+ if (S == nullptr) {
+ S = new double[min_mn]; // Singular values
+ owns_S = true;
}
- if (owns_S && S != nullptr) {
- delete[] S;
- S = nullptr;
+ if (VT == nullptr) {
+ VT = new double[min_mn * n]; // VT is min(m,n) x n
+ owns_VT = true;
}
- if (owns_VT && VT != nullptr) {
- delete[] VT;
- VT = nullptr;
- }
- 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;
+ 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;
- m = 0;
- n = 0;
- algo = Algorithm::GELSD;
- lambda = 0.0;
+ 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;
+ if (m != m_ || n != n_ || algo != algorithm) {
+ return false;
+ }
+ return true;
}
diff --git a/tests/test_ols.cpp b/tests/test_ols.cpp
index 8275369dc3fa029203504d580737b938503a4ee4..98a34d02998b5b81a2eb4bfa3527b7a8b4e1287b 100644
--- a/tests/test_ols.cpp
+++ b/tests/test_ols.cpp
@@ -18,8 +18,8 @@ TEST_CASE("Testing OLS") {
// 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
+ 2, 5, 8, // Second column
+ 3, 6, 10}); // Third column
Vector expected_weights = {1, 2, 3};
@@ -32,7 +32,51 @@ TEST_CASE("Testing OLS") {
vectors_are_close(weights, expected_weights, 1e-4);
}
- // Other sections remain unchanged...
+ 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
@@ -82,7 +126,7 @@ TEST_CASE("Testing OLS") {
// Recompute B_over
Vector B_over = A_over * expected_weights_over;
- Vector B_over2 = B_over;;
+ Vector B_over2 = B_over;
Vector weights(n);
@@ -121,12 +165,204 @@ TEST_CASE("Testing OLS") {
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);
+ 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;
+ };
+
+ // 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;
+ };
+
+ // 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);
+ double rcond = -1.0;
+ double lambda = 0.0005; // Regularization parameter
+
+ OLS::solve_with_svd(m, n, B, weights, svd, rcond, lambda);
+
+ // Since lambda introduces regularization, weights may not match exactly
+ // but for this test, we can still check that the solution is reasonable
+
+ Vector B_computed = A * weights;
+
+ // 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);
+ }
+ residual_norm = std::sqrt(residual_norm);
+
+ REQUIRE(residual_norm < 1e-3); // Adjust tolerance as needed
+ }
+
+ 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});
+
+ Vector B_ill = {1e-11, 1e-11, 1e-11};
+ Vector weights(3);
+
+ // 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("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);
+
+ OLS::solve(A_well, B_well, weights, -1.0, OLS::Algorithm::Cholesky);
+ vectors_are_close(weights, expected_weights, 1e-4);
+ }
+
// Rest of the tests...
}