EigValues.cpp 3.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110
  1. /**
  2. * @file EigValues.cpp
  3. * @brief EigValues Class
  4. * @author Michael Koch
  5. * @date 08/19/2008
  6. */
  7. #include <iostream>
  8. #include "EigValues.h"
  9. #define DEBUG_ARNOLDI
  10. using namespace NICE;
  11. using namespace std;
  12. void
  13. EVArnoldi::getEigenvalues ( const GenericMatrix & data, Vector & eigenvalues,
  14. Matrix & eigenvectors, uint k )
  15. {
  16. if ( data.rows () != data.cols () )
  17. {
  18. throw ( "EVArnoldi: matrix has to be quadratic" );
  19. }
  20. if ( k <= 0 )
  21. {
  22. throw ( "EVArnoldi: please use k>0." );
  23. }
  24. if ( verbose )
  25. cerr << "Initialize Matrices";
  26. uint n = data.cols ();
  27. NICE::Matrix rmatrix ( n, k, 0 ); //=eigenvectors
  28. NICE::Matrix qmatrix ( n, k, 0 ); //saves data =eigenvectors
  29. eigenvalues.resize ( k );
  30. NICE::Vector q ( n );
  31. NICE::Vector r ( n );
  32. if ( verbose )
  33. cerr << "Random Initialization" << endl;
  34. //random initialisation
  35. for ( uint i = 0; i < k; i++ )
  36. for ( uint j = 0; j < n; j++ )
  37. rmatrix ( j, i ) = drand48 ();
  38. // rmatrix ( j, i ) = 0.5;
  39. //TODO the random initialization might help, but it is bad for reproducibility :(
  40. //reduceddim
  41. double delta = 1.0;
  42. uint iteration = 0;
  43. while ( delta > mindelta && iteration < maxiterations )
  44. {
  45. //replace Vector rold by matrix rold to check for convergence of all eigenvectors
  46. //NICE::Vector rold ( rmatrix.getColumn ( k - 1 ) );
  47. NICE::Matrix rold(rmatrix);
  48. // meta-comment: i is an index for the iteration, j is an index for the basis
  49. // element (1 <= j <= k)
  50. for ( uint reduceddim = 0; reduceddim < k; reduceddim++ )
  51. {
  52. // -> get r^i_j from R matrix
  53. q = rmatrix.getColumn ( reduceddim );
  54. // q^i_j = r^i_j / ||r^i_j||
  55. q.normalizeL2 ();
  56. // -> store in Q matrix
  57. qmatrix.getColumnRef ( reduceddim ) = q;
  58. // this line copies a vector with external memory!
  59. // changing currentcol leads to a change in the R matrix!!
  60. Vector currentCol = rmatrix.getColumnRef ( reduceddim );
  61. // r^{i+1}_j = A * q^i_j ( r^i_j is overwritten by r^{i+1}_j )
  62. data.multiply ( currentCol, q );
  63. // for all j: r^{i+1}_j -= q^i_j * < q^i_j, r^{i+1}_j >
  64. for ( uint j = 0; j < reduceddim; j++ )
  65. rmatrix.getColumnRef ( reduceddim ) -=
  66. qmatrix.getColumn ( j ) *
  67. ( qmatrix.getColumn ( j ).
  68. scalarProduct ( rmatrix.getColumn ( reduceddim ) ) );
  69. }
  70. //convergence stuff (replaced by checking all eigenvectors instead of a single one
  71. //NICE::Vector diff = rold - rmatrix.getColumn ( k - 1 );
  72. //delta = diff.normL2 ();
  73. NICE::Vector tmpDiff;
  74. double norm_tmpDiff;
  75. delta = 0.0;
  76. for ( uint j = 0; j < k; j++ )
  77. {
  78. tmpDiff = rold.getColumn(j) - rmatrix.getColumn(j);
  79. norm_tmpDiff = tmpDiff.normL2();
  80. if (norm_tmpDiff > delta)
  81. delta = norm_tmpDiff;
  82. }
  83. iteration++;
  84. if ( verbose )
  85. cerr << "EVArnoldi: [" << iteration << "] delta=" << delta << endl;
  86. }
  87. eigenvectors = rmatrix;
  88. for ( uint i = 0; i < k; i++ )
  89. {
  90. NICE::Vector tmp;
  91. eigenvectors.getColumnRef ( i ).normalizeL2 ();
  92. data.multiply ( tmp, eigenvectors.getColumn ( i ) );
  93. eigenvalues[i] = tmp.scalarProduct ( eigenvectors.getColumn ( i ) );
  94. }
  95. }