Сопряженный градиент в собственных для эрмитовых матриц - PullRequest
Новая
       10

Сопряженный градиент в собственных для эрмитовых матриц

0 голосов
/ 05 мая 2018

в литературе метод сопряженных градиентов обычно представлен для реальных симметричных положительно определенных матриц. Однако в описании метода CG в библиотеке Eigen: https://eigen.tuxfamily.org/dox/group__IterativeLinearSolvers__Module.html можно найти утверждение: «Сопряженный градиент для самосопряженных (эрмитовых) матриц»

Это подразумевает, что он также должен работать для эрмитовых (сложных, а не чисто вещественных) матриц. Это тот случай?

Минимальный пример показывает, что на самом деле он не работает наивно с эрмитовыми матрицами. Есть ли хитрость, которую нужно знать, или это ошибка в описании?

В моем минимальном примере используются матрицы спина 3/2 Sx (вещественная симметрия) и Sy (комплексный эрмит), чьи собственные значения известны как -1,5, -0,5,0,5,1,5.

Результаты для реального симметричного случая хороши, но в сложном случае это приводит к NaN. 

#include <iostream> 
#include <complex>
#include <Eigen/Core>
#include <Eigen/IterativeLinearSolvers>

int main(int args, char **argv){

  Eigen::VectorXcd b=Eigen::VectorXcd::Ones(4);    
  Eigen::VectorXcd x;

  std::complex<double> i_unit(0,1);

//Hermitian matrix:
  Eigen::MatrixXcd A(4,4);
  A<<0,-i_unit*sqrt(3.)/2., 0 ,0,  \
     i_unit*sqrt(3.)/2., 0 ,-i_unit, 0,\
     0,i_unit,0,-i_unit*sqrt(3.)/2.,\
     0,0,i_unit*sqrt(3.)/2.,0;

//Real symmetric matrix:
  Eigen::MatrixXcd B(4,4);
  B<<0,sqrt(3.)/2., 0 ,0,  \
     sqrt(3.)/2., 0 ,1, 0,\
     0,1,0,sqrt(3.)/2.,\
     0,0,sqrt(3.)/2.,0;

  Eigen::ConjugateGradient< Eigen::MatrixXcd, Eigen::Lower|Eigen::Upper> cg;
  cg.compute(A);
  x = cg.solve(b);

  std::cout<<"Hermitian matrix:"<<std::endl;
  std::cout<<"A: "<<std::endl<<A<<std::endl;
  std::cout<<"b: "<<std::endl<<b<<std::endl;
  std::cout<<"x: "<<std::endl<<x<<std::endl;
  std::cout<<"(b-A*x).norm(): "<<(b-A*x).norm()<<std::endl;
  std::cout<<"cg.error(): "<<cg.error()<<std::endl;
  std::cout<<std::endl;

  cg.compute(B);
  x = cg.solve(b);
  std::cout<<"Real symmetric matrix:"<<std::endl;
  std::cout<<"B: "<<std::endl<<B<<std::endl;
  std::cout<<"b: "<<std::endl<<b<<std::endl;
  std::cout<<"x: "<<std::endl<<x<<std::endl;
  std::cout<<"(b-B*x).norm(): "<<(b-B*x).norm()<<std::endl;
  std::cout<<"cg.error(): "<<cg.error()<<std::endl;
  std::cout<<std::endl;
  return 0;
}

1 Ответ

0 голосов
/ 05 мая 2018

эрмитова недостаточно, оно также должно быть положительно определенным , что не соответствует вашему случаю, поскольку ваша матрица имеет как положительные, так и отрицательные собственные значения. В любом случае, CG скорее предназначен для обработки очень больших разреженных матриц, для матрицы 4x4 лучше использовать плотное разложение. В вашем случае, LDLT будет хорошо.

...