# [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:= [email protected]@{{0,1,0},{2,0,2},{0,1,0}}//FullForm
Out//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