[PYTHON] Ein Memorandum zur Verwendung von eigen3

Überblick

Wir verfolgen eine willkürliche Präzisionsdiagonalisierung (dichte Matrix, spärliche Matrix) und probieren verschiedene Dinge aus. Ich habe gestern eine Bibliothek namens eigen getroffen, sie heute installiert und verwendet, also habe ich sie als Protokoll aufgezeichnet.

Als willkürliche Präzisionsdiagonalisierung implementiert Mpack die Diagonalisierung der symmetrischen und Hermite-Klassen dichter Matrizen, aber es ist nicht möglich, grobe Matrizen effizient zu lösen, da spärliche Matrizen überhaupt nicht durch Lapack bereitgestellt werden. Es scheint, dass.

Nach vielen Recherchen scheint es, dass eine c ++ - Vorlagenbibliothek namens Eigen eine Diagonalisierung von dichten und groben Matrizen durchführen kann. Und wenn Sie die nicht unterstützte mpfr-Bibliothek verwenden, können Sie sie anscheinend mit beliebiger Genauigkeit diagonalisieren. Vergessen Sie also nicht, ein Protokoll mit eigen zu hinterlassen.

Der Hauptgrund für die Protokollierung ist, dass ich kein C ++ - Benutzer bin und es wahrscheinlich bald vergessen werde.

Umgebung ubuntu 14.04

install

Laden Sie das Neueste von Get it auf der Hauptseite Eigen herunter und entpacken Sie es.

Ich verstehe nicht, was eine Vorlagenbibliothek überhaupt ist, aber wenn ich Eigen Getting Start lese, muss die Vorlagenbibliothek installiert werden. Es scheint nicht zu existieren, und es scheint, dass es nur durch Ausführen des Programms im erweiterten Verzeichnis funktioniert.

Beispielprogramm in Erste Schritte in dem für die Testversion extrahierten Verzeichnis (das Verzeichnis mit Eigen und 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;
}

Ich habe versucht zu kompilieren

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

Weil es nicht kompiliert wurde

#include "Eigen/Core"

(Ich denke, es gab ein ähnliches Phänomen in C, aber ich habe es vergessen).

Eigens nicht unterstütztes Modul MPFRC ++ Support-Modul für weitere willkürliche Präzisionsberechnungen Stichprobe


#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;
}

Fügen Sie wie unten geändert hinzu

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

Und kompilieren

$ 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.

Aber es gibt einen Fehler. Da es keine Hilfe dafür gibt, kann das Beispielprogramm durch Kopieren von Eigen / und nicht unterstützt / unter / usr / local / include / ausgeführt werden.

MPRealSupport

MPRealSupport, da Sie beim Kompilieren die Bibliotheken gmp und mpfr verknüpfen müssen. Beispiele müssen wie folgt zusammengestellt werden

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

Vorbereitung für eine beliebige Präzisionsdiagonalisierung

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

Der Eigenwert von ist $ \ pm2,0 $. Wenn Sie mit doppelter Genauigkeit lösen, zum Beispiel numpy (Python)

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

Der Fehler breitet sich auf den 0-Eigenwert aus. Gleiches gilt, wenn Mathematica zur Diagonalisierung von maschinengenauen Matrizen verwendet wird.

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

Ist. Natürlich kann nMathematica mit beliebiger Genauigkeit diagonalisiert werden ($ \ le \ infty $) (aber ich mache mir Sorgen um die Geschwindigkeit). Programme, die mit eigen3 und mpfrc ++ diagonalisiert werden

#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;
}

Es scheint, dass Sie es so schreiben sollten, und das Ausführungsergebnis ist

$ 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)

Es ist ersichtlich, dass es mit hoher Genauigkeit ausgewertet wird. Wie in der Beschreibung

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

Dann ist das Ausführungsergebnis

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

Es ist nur so genau wie doppelt.

Interessanterweise stimmen das Ergebnis von Lapack (Numpy) mit 0 Eigenwert und das Ergebnis von Mathematica mit Maschinengenauigkeit (Double) im Genauigkeitsbereich überein. Wenn ich mich richtig erinnere, denke ich, dass Mathematica per Lapack implementiert wird, also wird es das gleiche Ergebnis sein.

Andererseits hat eigen3, soweit ersichtlich, unterschiedliche Werte, obwohl ich nicht weiß, wie ich die Anzahl der angezeigten Ziffern erhöhen kann. Liegt das an den unterschiedlichen Algorithmen?

Übrigens, wenn Sie mpmath (Python) verwenden, um mit der Standardgenauigkeit zu diagonalisieren (= 16 Stellen des falschen Teils)

>>> 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')],

Dies ist interessant, da auch die Ergebnisse unterschiedlich sind.

to do

Es ist einen halben Tag her, seit ich angefangen habe, C ++ zu verwenden. Ich weiß nicht, wie ich es benutzen soll.

Auch Geschwindigkeitsvergleich

(Zusatz)

Vergleichen Sie die Geschwindigkeiten bei Diagonalisierung mit eigen, mpack und mathematica. Da mpack eine begrenzte Klasse von Matrizen hat, die gelöst werden können, ist die folgende symmetrische Matrix

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}

Wird als dichte Matrix mit einer Genauigkeit von 80 Stellen diagonalisiert.

Immerhin ist Mpack früh dran. Da die CPU beim Ausführen des mpack-Programms 400 wird, kann sie durch blas parallelisiert werden. Eigen scheint in der Lage zu sein, Blas zu verwenden, aber was soll ich tun?

Die Legende dito (Anorldi 20%) ist das Eigensystem, eine diagonale Routine der Mathematik, die die Anorldi-Methode (Lanczos-Methode) spezifiziert und den um 20% größeren Eigenwert ermittelt.

Dies ist der Vergleich, wenn er aus dem kleineren berechnet wird Mit der Arnoldi-Methode und der Lanczos-Methode wird die Diagonalisierung von $ 10 ^ 5 \ times10 ^ 5 $ Realität.