[PYTHON] A memorandum of using eigen3

Overview

We pursue arbitrary precision diagonalization (dense matrix, sparse matrix) and try various things. I met a library called eigen yesterday, installed it today and used it, so I recorded it as a log.

As arbitrary precision diagonalization, Mpack implements diagonalization of symmetric and Hermite classes of dense matrices, but sparse matrices cannot be solved efficiently because lapack does not provide them in the first place. It seems.

After a lot of research, it seems that a c ++ template library called Eigen can diagonalize dense and coarse matrices. And if you use the unsupported mpfr library, it seems that you can diagonalize with arbitrary precision, so remember to leave a log using eigen.

The main reason for logging is that I'm not a c ++ user and I'm likely to forget it soon.

environment ubuntu 14.04

install

Download and extract latest from Get it on the main page Eigen.

I don't understand what a template library is in the first place, but if you read Eigen getting start, you need to make install the template library. It doesn't seem to exist, and if you run the program in the expanded directory, it seems to work by itself.

Sample program in getting start in the directory extracted for trial (the directory containing Eigen and 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;
}

I tried to compile

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

Because the compilation did not pass

#include "Eigen/Core"

(I think there was a similar phenomenon in C, but I forgot).

Eigen's unsupported module MPFRC ++ Support module for further arbitrary precision calculations Sample


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

Include modified as below

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

And compile

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

But it gives an error. Since there is no help for it, I made it possible to execute the sample program by sudo copying Eigen / and unsupported / under / usr / local / include /.

MPRealSupport

MPRealSupport because you need to link the gmp and mpfr libraries when compiling. Samples need to be compiled as follows

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

Preparation for arbitrary precision diagonalization

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

The eigenvalues of are $ \ pm2,0 $. If you solve it with double precision, for example numpy (python)

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

The error propagates to the 0 eigenvalues. The same is true when diagonalizing a machine-precision matrix in Mathematica.

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

Is. Of course, nMathematica can be diagonalized with arbitrary precision ($ \ le \ infinty $) (but I'm worried about speed). The program to diagonalize using eigen3 and 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;
}

It seems that you should write it like this, and the execution result is

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

It can be seen that it is evaluated with high accuracy. As in the description

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

Then the execution result is

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

It is only as accurate as double.

Interestingly, the 0-eigenvalue lapack (numpy) results and the machine-precision (double) Mathematica results agree in the precision range. If I remember correctly, I think Mathematica is implemented by lapack, so it will be the same result.

On the other hand, eigen3 has different values as far as it can be seen, although I do not know how to increase the number of displayed digits. Is this due to the different algorithms used?

By the way, if you diagonalize with default precision (= 16 digits of mantissa) with mpmath (python)

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

This is interesting because the results are also different.

to do

Is MPRealSupport only real support by name? Is complex numbers also supported? It's been half a day since I started using c ++ in the first place. I don't know how to use it.

Also speed comparison

(Addition)

Compare the velocities when diagonalized with eigen, mpack, and mathematica. Since mpack has a limited class of matrices that can be solved, the following symmetric 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}

Is diagonalized as a dense matrix with a precision of 80 digits.

diag-time.png

After all Mpack is early. Since the cpu becomes 400 while executing the mpack program, it may be parallelized by blas. Eigen seems to be able to use blas, but what should I do?

The legend ditto (Anorldi 20%) is the diagonalization routine of mathematica, Eigensystem, which specifies the Anorldi method (Lanczos method) and finds the larger eigenvalue by 20%.

This is the comparison when calculated from the smaller one diag-time-arnoldi.png With the Arnoldi method and Lanczos method, the diagonalization of $ 10 ^ 5 \ times10 ^ 5 $ becomes realistic.

Recommended Posts

A memorandum of using eigen3
A memorandum of using Python's input function
A memorandum of kernel compilation
A small memorandum of openpyxl
[Python] A memorandum of beautiful soup4
A memorandum of files under conf.d
A memorandum of closure survey contents
Getting a combination of elements using itertools
Memorandum of sed
A memorandum of speed of arbitrary degree diagonalization
Impressions of using Flask for a month
A memorandum of understanding about django's QueryDict
A memorandum of python string deletion process
A memorandum of trouble when formatting data
A memorandum of calling Python from Common Lisp
A memorandum of studying and implementing deep learning
A memorandum of extraction by python bs4 request
[Linux command] A memorandum of frequently used commands
Memorandum of fastText (editing)
memorandum of vi command
A memorandum about matplotlib
A memorandum about Nan.
elasticsearch_dsl Memorandum of Understanding
Example of using lambda
A memorandum of method often used in machine learning using scikit-learn (for beginners)
[Python] Implementation of clustering using a mixed Gaussian model
A memorandum about the warning of the pylint output result
Cut a part of the string using a Python slice
Read a large amount of securities reports using COTOHA
I tried using Python (3) instead of a scientific calculator
Instantly create a diagram of 2D data using python's matplotlib
Implementation of VGG16 using Keras created without using a trained model
Avoiding the pitfalls of using a Mac (for Linux users?)
Implementation of TF-IDF using gensim
Reuse the behavior of the @property method by using a descriptor [16/100]
Time measurement using a clock
python: Basics of using scikit-learn ①
Pepper Tutorial (5): Using a Tablet
A memorandum about correlation [Python]
[Introduction to AWS] A memorandum of building a web server on AWS
Using a printer with Debian 10
A memo of writing a basic function in Python using recursion
A memorandum about Python mock
A memorandum regarding γ conversion
A brief summary of Linux
The story of creating a database using the Google Analytics API
[End of 2020] A memo to start using AWS CLI (Version 2)
A memorandum of understanding for the Python package management tool ez_setup
I tried to get a database of horse racing using Pandas
I tried to make a regular expression of "amount" using Python
I tried to make a regular expression of "time" using Python
A memorandum until using mecab on a machine that cannot use sudo
I tried to make a regular expression of "date" using Python
Convert a large number of PDF files to text files using pdfminer
I tried to get a list of AMI Names using Boto3
How to save only a part of a long video using OpenCV
Evaluate the performance of a simple regression model using LeaveOneOut cross-validation
Output search results of posts to a file using Mattermost API
A memorandum of scraping & machine learning [development technique] by Python (Chapter 4)
A memorandum of scraping & machine learning [development technique] by Python (Chapter 5)
Try using Elasticsearch as the foundation of a question answering system