[PYTHON] Un mémorandum d'utilisation de eigen3

Aperçu

Nous poursuivons une diagonalisation de précision arbitraire (matrice dense, matrice clairsemée) et essayons diverses choses. J'ai rencontré une bibliothèque appelée eigen hier, je l'ai installée aujourd'hui et je l'ai utilisée, alors je l'ai enregistrée sous forme de journal.

En tant que diagonalisation de précision arbitraire, Mpack implémente la diagonalisation des classes symétrique et Hermite des matrices denses, mais il n'est pas possible de résoudre efficacement des matrices grossières car les matrices clairsemées ne sont pas fournies par lapack en premier lieu. Il semble.

Après de nombreuses recherches, il semble qu'une bibliothèque de modèles c ++ appelée Eigen puisse effectuer la diagonalisation de matrices denses et grossières. Et si vous utilisez la bibliothèque mpfr non prise en charge, il semble que vous puissiez la diagonaliser avec une précision arbitraire, alors n'oubliez pas de laisser un journal en utilisant eigen.

La principale raison de la journalisation est que je ne suis pas un utilisateur C ++ et que je vais probablement l'oublier bientôt.

environnement ubuntu 14.04

install

Téléchargez et décompressez la dernière version de Get it sur la page principale Eigen.

Je ne comprends pas ce qu'est une bibliothèque de modèles en premier lieu, mais quand je lis Eigen getting start, la bibliothèque de modèles doit être installée. Il ne semble pas exister, et si vous exécutez le programme dans le répertoire développé, il semble fonctionner tout seul.

Exemple de programme dans getting start dans le répertoire extrait pour l'essai (le répertoire où se trouvent Eigen et README.md)

#include <iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
int main()
{
  MatrixXd m(2,2);
  m(0,0) = 3;
  m(1,0) = 2.5;
  m(0,1) = -1;
  m(1,1) = m(1,0) + m(0,1);
  std::cout << m << std::endl;
}

J'ai essayé de compiler

$ g++ test.cpp　-o a.out
test.cpp:3:23: fatal error: Eigen/Dense: No such file or directory
compilation terminated.

Parce qu'il n'a pas compilé

#include "Eigen/Core"

(Je pense qu'il y avait un phénomène similaire en C, mais j'ai oublié).

Module non pris en charge d'Eigen Module de support MPFRC ++ pour d'autres calculs de précision arbitraire Échantillon


#include <iostream>
#include <Eigen/MPRealSupport>
#include <Eigen/LU>
using namespace mpfr;
using namespace Eigen;
int main()
{
  // set precision to 256 bits (double has only 53 bits)
  mpreal::set_default_prec(256);
  // Declare matrix and vector types with multi-precision scalar type
  typedef Matrix<mpreal,Dynamic,Dynamic>  MatrixXmp;
  typedef Matrix<mpreal,Dynamic,1>        VectorXmp;
  MatrixXmp A = MatrixXmp::Random(100,100);
  VectorXmp b = VectorXmp::Random(100);
  // Solve Ax=b using LU
  VectorXmp x = A.lu().solve(b);
  std::cout << "relative error: " << (A*x - b).norm() / b.norm() << std::endl;
  return 0;
}

Inclure modifié comme ci-dessous

#include "unsupported/Eigen/MPRealSupport"
#include "Eigen/LU"

Et compilez

$ g++ test.cpp -o a.out -lgmp -lmpfr
In file included from test.cpp:2:0:
unsupported/Eigen/MPRealSupport:15:22: fatal error: Eigen/Core: No such file or directory
compilation terminated.

Mais cela donne une erreur. Puisqu'il n'y a aucune aide pour cela, j'ai rendu possible l'exécution de l'exemple de programme en copiant Eigen / et unsupported / sous / usr / local / include /.

MPRealSupport

MPRealSupport car vous devez lier les bibliothèques gmp et mpfr lors de la compilation. Les échantillons doivent être compilés comme suit

g++ test.cpp -o a -lgmp -lmpfr

Préparation à une diagonalisation de précision arbitraire

queue $ A = \begin{pmatrix} 0 & 1 & 0\\\\ 2 & 0 & 2\\\\ 0 & 1 & 0 \end{pmatrix} $

La valeur propre de est $ \ pm2,0 $. Si vous résolvez avec une double précision, par exemple numpy (python)

>>> import numpy as np
>>> np.set_printoptions(precision=18)
array([ -2.000000000000000000e+00,  -1.732555071631836086e-16,
         1.999999999999999556e+00])

L'erreur se propage à la valeur propre 0. Il en va de même lorsque Mathematica est utilisé pour diagonaliser des matrices de précision machine.

In[15]:= Eigenvalues@N@{{0,1,0},{2,0,2},{0,1,0}}//FullForm
Out[15]//FullForm= List[-2.`,1.9999999999999996`,-1.732555071631836`*^-16]

Est. Bien sûr, nMathematica peut être diagonalisé avec une précision arbitraire ($ \ le \ infty $) (mais je m'inquiète pour la vitesse). Programmes diagonalisés en utilisant eigen3 et mpfrc ++

#include <iostream>
#include <Eigen/Eigenvalues>
#include <unsupported/Eigen/MPRealSupport>
#include <Eigen/Dense>
using namespace Eigen;

using namespace mpfr;

int main()
{
  // set precision to 256 bits (double has only 53 bits)
  mpreal::set_default_prec(53);
  // Declare matrix and vector types with multi-precision scalar type
  typedef Matrix<mpreal,Dynamic,Dynamic>  MatrixXmp;
  MatrixXmp A(3,3); 

  A(0,0) = 0;
  A(0,1) = 1;
  A(0,2) = 0;

  A(1,0) = 2;
  A(1,1) = 0;
  A(1,2) = 2;

  A(2,0) = 0;
  A(2,1) = 1;
  A(2,2) = 0;
  
  std::cout << A << std::endl;
    
  EigenSolver<MatrixXmp> es(A);
  std::cout << es.eigenvalues() << std::endl;

  return 0;
}

Il semble que vous devriez l'écrire comme ceci, et le résultat de l'exécution est

$ g++ eigenvalues-mpfrc-min.cpp -o a -lgmp -lmpfr 
$ ./a 
0 1 0
2 0 2
0 1 0
          (-2,0)
           (2,0)
(-3.28269e-77,0)

On voit qu'il est évalué avec une grande précision. Comme dans la description

  // set precision to 256 bits (double has only 53 bits)
  mpreal::set_default_prec(53);

Alors le résultat de l'exécution est

$ ./a 
0 1 0
2 0 2
0 1 0
         (-2,0)
          (2,0)
(1.77435e-16,0)

C'est seulement aussi précis que le double.

Fait intéressant, le résultat de lapack (numpy) avec une valeur propre 0 et le résultat de Mathematica avec une précision machine (double) sont en accord dans la plage de précision. Si je me souviens bien, je pense que Mathematica est implémenté par lapack, donc ce sera le même résultat.

Par contre, eigen3 a des valeurs différentes pour autant qu'il soit visible, même si je ne sais pas comment augmenter le nombre de chiffres affichés. Est-ce dû aux différents algorithmes utilisés?

À propos, si vous utilisez mpmath (python) pour diagonaliser avec la précision par défaut (= 16 chiffres de la partie incorrecte)

>>> import mpmath
>>> mpmath.eig(mpmath.matrix([[0,1,0],[2,0,2],[0,1,0]]))
([mpf('-1.9999999999999993'), mpf('1.6358461123585456e-16'), mpf('2.0')],

Ceci est intéressant car les résultats sont également différents.

to do

MPRealSupport est-il uniquement un support réel du nom? Les nombres complexes sont-ils également pris en charge? Cela fait une demi-journée depuis que j'ai commencé à utiliser C ++. Je ne sais pas comment m'en servir.

Aussi la comparaison de vitesse

(Une addition)

Comparez les vitesses en diagonale avec eigen, mpack et mathematica. Puisque mpack a une classe limitée de matrices qui peuvent être résolues, la matrice symétrique suivante

M= \begin{pmatrix} 1&1&0&\cdots &0\\\\ 1&1&1&\cdots &0\\\\ 0&1&1&\cdots &0\\\\ &\vdots&&\ddots&\\\\ 0&0&0&\cdots &1\\\\ \end{pmatrix}

Est diagonalisée comme une matrice dense avec une précision de 80 chiffres.

Après tout, Mpack est en avance. Puisque cpu devient 400 lors de l'exécution du programme mpack, il peut être parallélisé par blas. Eigen semble être capable d'utiliser des blas, mais que dois-je faire?

La légende idem (Anorldi 20%) est le système Eigen, qui est une routine diagonale de Mathematica, qui spécifie la méthode Anorldi (méthode Lanczos) et trouve la plus grande valeur propre de 20%.

C'est la comparaison calculée à partir du plus petit Avec la méthode Arnoldi et la méthode de Lanczos, la diagonalisation de $ 10 ^ 5 \ times10 ^ 5 $ devient une réalité.