viennacl / viennacl-dev

Developer repository for ViennaCL. Visit http://viennacl.sourceforge.net/ for the latest releases.
Other
281 stars 89 forks source link

Tuning Lanczos to find all eigenvalues of a dense Hermitian matrix #280

Open TysonRayJones opened 5 years ago

TysonRayJones commented 5 years ago

I'm encoding an n-by-n Hermitian matrix

[ a + b i   ... ]
[ ...           ]

into a 2n-by-2n ViennaCL matrix

[ a  -b  ... ]
[ b   a  ... ]
[ ...        ]

which mathematically, should have the same (though each twice-degenerate) real eigenvalues.

I'm having trouble tuning Lanczos method to accurately compute these eigenvalues, using my Quadro P6000 (compute capability 6.1).

E.g. this 4-by-4 Hermitian matrix (in Mathematica):

{{0.332082 + 0. I, 0.669948 + 0.405974 I, 0.677182 + 0.194996 I, 
  0.13695 + 0.148397 I}, {0.669948 - 0.405974 I, 0.0106002 + 0. I, 
  0.330906 + 0.132972 I, 
  0.747063 + 0.0849097 I}, {0.677182 - 0.194996 I, 
  0.330906 - 0.132972 I, 0.690913 + 0. I, 
  0.352845 + 0.00661738 I}, {0.13695 - 0.148397 I, 
  0.747063 - 0.0849097 I, 0.352845 - 0.00661738 I, 0.880548 + 0. I}}

with eigenvalues {-0.835489, 0.0414574, 0.671312, 2.03686}, becomes this ViennaCL matrix:

viennacl::matrix<double> m(8,8);
m(0,0)=0.3320820211; m(0,1)=0.; m(0,2)=0.6699477692; m(0,3)=-0.405973848; m(0,4)=0.677182054; m(0,5)=-0.1949962063; m(0,6)=0.1369502543; m(0,7)=-0.1483970803;
m(1,0)=0.; m(1,1)=0.3320820211; m(1,2)=0.405973848; m(1,3)=0.6699477692; m(1,4)=0.1949962063; m(1,5)=0.677182054; m(1,6)=0.1483970803; m(1,7)=0.1369502543;
m(2,0)=0.6699477692; m(2,1)=0.405973848; m(2,2)=0.01060019334; m(2,3)=0.; m(2,4)=0.3309064604; m(2,5)=-0.1329722856; m(2,6)=0.7470626654; m(2,7)=-0.08490966432;
m(3,0)=-0.405973848; m(3,1)=0.6699477692; m(3,2)=0.; m(3,3)=0.01060019334; m(3,4)=0.1329722856; m(3,5)=0.3309064604; m(3,6)=0.08490966432; m(3,7)=0.7470626654;
m(4,0)=0.677182054; m(4,1)=0.1949962063; m(4,2)=0.3309064604; m(4,3)=0.1329722856; m(4,4)=0.6909132009; m(4,5)=0.; m(4,6)=0.3528451954; m(4,7)=-0.006617375307;
m(5,0)=-0.1949962063; m(5,1)=0.677182054; m(5,2)=-0.1329722856; m(5,3)=0.3309064604; m(5,4)=0.; m(5,5)=0.6909132009; m(5,6)=0.006617375307; m(5,7)=0.3528451954;
m(6,0)=0.1369502543; m(6,1)=0.1483970803; m(6,2)=0.7470626654; m(6,3)=0.08490966432; m(6,4)=0.3528451954; m(6,5)=0.006617375307; m(6,6)=0.880548364; m(6,7)=0.;
m(7,0)=-0.1483970803; m(7,1)=0.1369502543; m(7,2)=-0.08490966432; m(7,3)=0.7470626654; m(7,4)=-0.006617375307; m(7,5)=0.3528451954; m(7,6)=0.; m(7,7)=0.880548364;

which I attempt to eigendecompose with

viennacl::linalg::lanczos_tag ltag(0.9,  8, 1, 2*8);
std::vector<double> lanczos_eigenvalues = viennacl::linalg::eig(m, ltag);

for (std::size_t i = 0; i < lanczos_eigenvalues.size(); i++)
        std::cout << "Approx. eigenvalue " << std::setprecision(7) << lanczos_eigenvalues[i] << "\n";

The eigenvalues reported by ViennaCL not only vary drastically in accuracy based on the parameters to ltag (which I'm unsure how to optimise), but vary between runs with the same parameters! Here are some outputs without changing the code:

# run 1
Approx. eigenvalue 2.036864
Approx. eigenvalue 2.034437
Approx. eigenvalue 0.6713113
Approx. eigenvalue 0.6667219
Approx. eigenvalue 0.04680611
Approx. eigenvalue 0.04145725
Approx. eigenvalue -0.8302918
Approx. eigenvalue -0.8354891
# run 2
Approx. eigenvalue 2.036864
Approx. eigenvalue 2.02396
Approx. eigenvalue 0.7262711
Approx. eigenvalue 0.6713118
Approx. eigenvalue 0.05332146
Approx. eigenvalue 0.04145729
Approx. eigenvalue -0.77557
Approx. eigenvalue -0.835489
# run 3
Approx. eigenvalue 2.036864
Approx. eigenvalue 2.036864
Approx. eigenvalue 0.6713121
Approx. eigenvalue 0.6713115
Approx. eigenvalue 0.04145736
Approx. eigenvalue 0.04145736
Approx. eigenvalue -0.8354883
Approx. eigenvalue -0.8354883

As you can tell just by comparing the degenerate pairs of eigenvalues, run 1 is inaccurate, run 2 is wildly inaccurate, and run 3 is excellent.

I need to configure ltag so that, no matter the runtime, the eigenvalues are computed very accurately every time. My end-goal is to accurately compute the sum of the absolute value of all the eigenvalues of Hermitian matrices as large as ~30k-by-30k!

Is there any elaboration on the parameters only briefly mentioned here, and how they impact the runtime and final precision?

TysonRayJones commented 5 years ago

Is there any doc on tuning Lanczos in ViennaCL I could grab a quick link to? Otherwise I'll need to switch to another library :(

karlrupp commented 4 years ago

@TysonRayJones Unfortunately there is no good documentation on Lanczos in ViennaCL. This is contributed code, so I don't know all the implementation details. Sorry!