dashSVD.h

Go to the documentation of this file.

 
// /************************************************************
//  * @file        dashSVD.h
//  * @author      Gabriel Hoffman
//  * @email       gabriel.hoffman@mssm.edu
//  * @brief       dashSVD One algorithm
//  * Copyright (C) 2026 Gabriel Hoffman
//  *************************************************/
 
 
// #include "GenomicDataStream.h"
 
// using namespace Rcpp;
// using namespace arma;
 
// struct SVDResult {
//   mat U;
//   vec S;
//   mat V;
// };
 
 
// struct PCA {
 
//   PCA( mat u):
//     u(u) {}
 
//   PCA( const mat &u, const mat &v, const vec &d, const vector<string> featureIds):
//     u(u), v(v), d(d), featureIds(featureIds) {}
 
//   mat u, v;
//   vec d;
//   const vector<string> featureIds;
// };
 
 
// SVDResult eigSVD_internal(const mat& A, double tol = 1e-10) {
//   mat C = A.t() * A;
 
//   vec eigval;
//   mat V;
 
//   eig_sym(eigval, V, C);
 
//   eigval = reverse(eigval);
//   V = fliplr(V);
 
//   vec sigma = sqrt(clamp(eigval, 0.0, datum::inf));
 
//   vec sigma_inv = zeros<vec>(sigma.n_elem);
//   uvec keep = find(sigma > tol);
//   sigma_inv.elem(keep) = 1.0 / sigma.elem(keep);
 
//   mat U = A * V;
//   U.each_row() %= sigma_inv.t();
 
//   return {U, sigma, V};
// }
 
 
// inline mat mat_prod_AtAQ(const mat& A, const mat& Q) {
//   return A.t() * (A * Q);
// }
 
// inline mat crossprod(const mat& A, const mat& Q) {
//   return A.t() * Q;
// }
 
 
// // [[Rcpp::export]]
// List dashSVD(
//     const mat& A,
//     int k,
//     int p = 3,
//     int s = 20,
//     int rand = 1,
//     bool verbose = false,
//     double tol = 1e-10
// ) {
//   const uword M = A.n_rows;
//   const uword N = A.n_cols;
//   const uword L = static_cast<uword>(k + s);
 
//   if (k <= 0) {
//     stop("k must be positive");
//   }
 
//   if (L > N) {
//     stop("k + s must be <= ncol(A)");
//   }
 
//   mat Omg;
 
//   if (rand == 1) {
//     Omg = randn<mat>(M, L);
//   } else {
//     Omg = randu<mat>(M, L);
//   }
 
//   mat Q = crossprod(A.t(), Omg);
 
//   SVDResult e = eigSVD_internal(Q, tol);
//   Q = e.U;
 
//   double alpha = 0.0;
 
//   for (int i = 1; i <= p; ++i) {
//     if (verbose) {
//       Rcpp::Rcout << "power iteration " << i
//                   << ", alpha=" << alpha << std::endl;
//     }
 
//     mat B = mat_prod_AtAQ(A, Q) - alpha * Q;
 
//     e = eigSVD_internal(B, tol);
//     Q = e.U;
 
//     if (e.S(L - 1) > alpha) {
//       alpha = (alpha + e.S(L - 1)) / 2.0;
//     }
//   }
 
//   e = eigSVD_internal(A * Q, tol);
 
//   mat U = e.U;
//   vec S = e.S;
//   mat V = Q * e.V;
 
//   const uword kk = static_cast<uword>(k);
 
//   return PCA(U.cols(0, kk - 1),
//             V.cols(0, kk - 1),
//             S.head(kk))
// }