Vespucci  1.0.0
svds.cpp
Go to the documentation of this file.
3 
24 bool Vespucci::Math::DimensionReduction::svds(const arma::mat &X, arma::uword k, arma::mat &U, arma::vec &s, arma::mat &V)
25 {
26  if (X.n_cols < 2){
27  std::cerr << "svds: X must be 2D arma::matrix" << std::endl;
28  return false;
29  }
30 
31  arma::uword m = X.n_rows;
32  arma::uword n = X.n_cols;
33  arma::uword p = (m < n ? m : n);
34  //arma::uword p = Vespucci::Math::min(m, n); //used below to establish tolerances.
35  //arma::uword q = Vespucci::Math::max(m, n);
36 
37  if (k > p){
38  if (k > m)
39  std::cerr << "svds: value of k " << k << "is greater than number of rows " << m << std::endl;
40  else
41  std::cerr << "svds: value of k" << k << "is greater than number of columns " << n << std::endl;
42  }
43 
44  k = (p < k ? p : k);
45  //Vespucci::Math::min(p, k); //make sure value of k is less than smallest dimension
46 
47  arma::sp_mat B(m+n, m+n);
48  B.submat(arma::span(0, m-1), arma::span(m, m+n-1)) = X; // top right corner
49  B.submat(arma::span(m, m+n-1), arma::span(0, m-1)) = X.t(); //bottom left corner
50  //If, for some reason, a arma::matrix of zeroes is input...
51  if (B.n_nonzero == 0){
52  U = arma::eye(m, k);
53  s = arma::zeros(k, k);
54  V = arma::eye(n, k);
55  std::cerr << "svds: input arma::matrix has no non-zero elements" << std::endl;
56  return false;
57  }
58 
59  arma::vec eigval;
60  arma::mat eigvec;
61 
62  //because we're using sigma=0, results of eigs will be centered around 0
63  //we throw out the negative ones, then add in 0 eigenvalues if we have to
64  //trying some weird stuff to figure out why arpack call fails
65  bool eigs_success = eigs_sym(eigval, eigvec, B, 2*k);
66 
67  if (!eigs_success){
68  std::cerr << "svds: arma::eigs_sym did not converge!" << std::endl;
69  }
70 
71  //Manipulate the results
72 
73  //sort eigvec from largest to smallest:
74  eigvec = eigvec.cols(sort_index(eigval, "descend"));
75 
76  //the negative eigenvalues are artifacts of how the matrix is constructed
77  //it is possible that there are two 0 eigenvalues. Only the first is taken
78  //The octave and matlab svds routines are a little more precise with this...
79  //It is highly unlikely that a 0 eigenvalue will occur with the types of data
80  //processed by Vespucci
81  arma::uvec return_indices = find(eigval >= 0, k, "first");
82 
83  U = sqrt(2) * (eigvec.rows(0, m-1));
84  U = U.cols(return_indices);
85 
86  s = eigval.rows(return_indices);
87 
88  V = sqrt(2) * eigvec.rows(m, m+n-1);
89  V = V.cols(return_indices);
90 
91  arma::uvec indices = sort_index(s, "descend");
92  s = sort(s, "descend");
93 
94  U = U.cols(indices);
95  V = V.cols(indices);
96  return eigs_success;
97 }
98 
VESPUCCI_EXPORT bool svds(const arma::mat &X, arma::uword k, arma::mat &U, arma::vec &s, arma::mat &V)
Vespucci::Math::DimensionReduction::svds Finds a few largest singular values of the arma::matrix X...
Definition: svds.cpp:24