52 throw std::runtime_error(
53 "Cannot calculate reduce_to_tridiagonal_matrix() with an empty matrix." 57 throw std::runtime_error(
58 "Expecting symmetrical matrix for reduce_to_tridiagonal_matrix()" 67 for(
size_t i = data.
cols()-1; i >= 1; --i ) {
68 size_t const l = i - 1;
74 for(
size_t k = 0; k <= l; ++k ) {
75 scale += std::fabs(data( i, k ));
81 for(
size_t k = 0; k <= l; ++k ) {
82 data( i, k ) /=
scale;
83 h += data( i, k ) * data( i, k );
86 double f = data( i, l );
87 double g = ( f > 0 ? -sqrt(h) : sqrt(h) );
93 for(
size_t j = 0; j <= l; ++j ) {
94 data( j, i ) = data( i, j )/h;
97 for(
size_t k = 0; k <= j; ++k ) {
98 g += data( j, k ) * data( i, k );
100 for(
size_t k = j+1; k <= l; ++k ) {
101 g += data( k, j ) * data( i, k );
108 double hh = f / (h + h);
109 for(
size_t j = 0; j <= l; ++j ) {
112 for(
size_t k = 0; k <= j; ++k ) {
113 data( j, k ) -= (f * tri.
intermediates[k] + g * data( i, k ));
127 for(
size_t i = 0; i < data.
cols(); ++i ) {
129 size_t const l = i - 1;
130 for(
size_t j = 0; j <= l; ++j ) {
132 for(
size_t k = 0; k <= l; ++k ) {
133 g += data( i, k ) * data( k, j );
135 for(
size_t k = 0; k <= l; ++k ) {
136 data( k, j ) -= g * data( k, i );
145 size_t const l = i - 1;
146 for(
size_t j = 0; j <= l; ++j ) {
163 size_t max_iterations
166 throw std::runtime_error(
167 "Cannot calculate tridiagonal_ql_algorithm() with an empty matrix." 171 throw std::runtime_error(
172 "Expecting symmetrical matrix for tridiagonal_ql_algorithm()" 176 throw std::runtime_error(
177 "Expecting TridiagonalDecompositionData vectors of the same size" 178 " as the data matrix in tridiagonal_ql_algorithm()" 185 size_t const n = data.
rows();
187 for(
size_t i = 1; i < n; ++i ) {
193 for(
size_t l = 0; l < n; ++l ) {
198 for( m = l; m < n-1; ++m ) {
199 double const dd = std::fabs( d[m]) + std::fabs(d[m+1] );
200 if( std::fabs(e[m]) + dd == dd ) {
206 if( max_iterations > 0 && iter == max_iterations ) {
207 throw std::runtime_error(
"No convergence in tridiagonal_ql_algorithm()." );
215 double g = (d[l+1] - d[l]) / (2.0 * e[l]);
216 double r = sqrt((g * g) + 1.0);
217 auto sign = ( (g) < 0 ? -std::fabs(r) : std::fabs(r) );
218 g = d[m] - d[l] + e[l] / (g + sign);
220 for(
long i = static_cast<long>( m ) - 1; i >=
static_cast<long>( l ); --i ) {
224 if( std::fabs(f) >= std::fabs(g) ) {
226 r = sqrt((c * c) + 1.0);
231 r = sqrt((s * s) + 1.0);
237 r = (d[i] - g) * s + 2.0 * c * b;
241 for(
size_t k = 0; k < n; ++k ) {
243 data( k, i+1 ) = s * data( k, i ) + c * f;
244 data( k, i ) = c * data( k, i ) - s * f;
275 auto standardized_data = data;
282 for(
auto& elem : symmat ) {
283 elem /=
static_cast<double>( standardized_data.rows() );
290 if( components == 0 ) {
291 components = standardized_data.cols();
293 if( components > standardized_data.cols() ) {
294 throw std::runtime_error(
295 "Cannot calculate more PCA components than the original data has columns." 304 assert( tri.eigenvalues.size() == standardized_data.cols() );
305 assert( tri.intermediates.size() == standardized_data.cols() );
306 assert( symmat.rows() == standardized_data.cols() );
307 assert( symmat.cols() == standardized_data.cols() );
311 tri.eigenvalues.begin(),
312 tri.eigenvalues.end(),
313 std::greater<double>()
315 result.
eigenvalues = std::vector<double>( components );
317 for(
size_t c = 0; c < components; ++c ) {
318 result.
eigenvalues[c] = tri.eigenvalues[ sorted_indices[c] ];
320 for(
size_t r = 0; r < symmat.rows(); ++r ) {
321 result.
eigenvectors( r, c ) = symmat( r, sorted_indices[c] );
Matrix< T > matrix_multiplication(Matrix< A > const &a, Matrix< B > const &b)
Calculate the product of two Matrices.
Provides some valuable algorithms that are not part of the C++ 11 STL.
std::vector< double > eigenvalues
Matrix< double > sums_of_squares_and_cross_products_matrix(Matrix< double > const &data)
Calculate the Sums of Squares and Cross Products Matrix (SSCP Matrix).
std::vector< double > intermediates
Standardize the mean, but not the variance of the data before performing the PCA. ...
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
Helper stucture that collects the output of principal_component_analysis().
std::vector< size_t > sort_indices(RandomAccessIterator first, RandomAccessIterator last, Comparator comparator)
Get the indices to the sorted order of the given range.
PcaData principal_component_analysis(Matrix< double > const &data, size_t components, PcaStandardization standardization)
Perfom a Principal Component Analysis on a given data Matrix.
Matrix< double > eigenvectors
TridiagonalDecompositionData reduce_to_tridiagonal_matrix(Matrix< double > &data)
Triangular decomposition of a symmetric matrix.
void scale(Histogram &h, double factor)
Matrix< double > projection
PcaStandardization
Setting for principal_component_analysis() to determine which form of standardization of the data to ...
void tridiagonal_ql_algorithm(Matrix< double > &data, TridiagonalDecompositionData &tri, size_t max_iterations)
Reduce a symmetric matrix to a symmetric tridiagonal matrix.
std::vector< double > eigenvalues
Standardize the mean and variance of the data before performing the PCA.
Helper structure used for the eigenvalue decomposition in reduce_to_tridiagonal_matrix() and tridiago...
std::vector< MeanStddevPair > standardize_cols(Matrix< double > &data, bool scale_means, bool scale_std)
Standardize the columns of a Matrix by subtracting the mean and scaling to unit variance.