Add expm operation #32
Comments
|
@masa16 I have implemented the expm method with the basic algorithm using the Padé approximation. Please review my pull request above. |
|
@seoanezonjic Hi. The expm method has been added to the master branch. Please install and confirm that. $ gem install specific_install
$ gem specific_install https://github.com/ruby-numo/numo-linalg.git masterNumo::Linalg.expm(a)
# The ord parameter value affects the approximation accuracy.
Numo::Linalg.expm(a, 20) |
|
Hi @yoshoku Two different chunks of code give me this error: diagonal_matrix = matrix.sum(1).diag
matrix_L = diagonal_matrix - matrix
beta = 0.02
beta_product = matrix_L * -beta
matrix_result = Numo::Linalg.expm(beta_product, 12)diagonal_matrix = matrix.sum(1).diag
matrix_L = diagonal_matrix - matrix
beta=0.04
id_mat = Numo::DFloat.eye(dimension_elements)
m_matrix = (id_mat * dimension_elements - diagonal_matrix + matrix ) * (beta/dimension_elements)
matrix_result = Numo::Linalg.expm(m_matrix, 12)matrix has 10 000 x 10 000 random float elements ranged between 0 and 100 but it's a symmetric matrix. What is the problem tha causes this seg fault? |
|
@seoanezonjic I think that the segmentation fault error occurs depending on the execution environment. That error did not occur in my environment (OS: macOS 10.15.1, CPU: Dual-Core Intel Core m3, Memory: 8 GB). Perhaps the cause of the error is in the BLAS/LAPACK library. I setup the environment using OpenBLAS for the BLAS/LAPACK library. $ brew install openblas
$ rbenv install 2.4.1
$ rbenv global 2.4.1
$ gem install specific_install
$ gem specific_install https://github.com/ruby-numo/numo-linalg.git masterThen, I executed the following code. require 'numo/linalg/autoloader'
puts 'generate symmetric matrix'
matrix = Numo::DFloat.new(10000, 10000).rand
matrix = 0.5 * (matrix + matrix.transpose)
puts 'calculate laplacian matrix'
diagonal_matrix = matrix.sum(1).diag
matrix_L = diagonal_matrix - matrix
beta = 0.02
beta_product = matrix_L * -beta
puts 'calculate matrix exponential'
matrix_result = Numo::Linalg.expm(beta_product, 12)$ ruby foo.rb
generate symmetric matrix
calculate laplacian matrix
calculate matrix exponential
$ |
|
Hi @yoshoku |
|
Hi @yoshoku |
|
@seoanezonjic I implemented the expm method with reference to the book [1]. How to find the order The The If the matrix norm and the acceptable approximation error are clear, you may be able to find def factorial(n)
n.zero? ? 1 : n * factorial(n-1)
end
def eps(q)
(2**(3 - 2 * q)) * (factorial(q)**2).fdiv(factorial(2 * q) * factorial(2 * q + 1))
end
puts eps(14) # 8.402086919997803e-47
puts eps(16) # 3.5687151066343315e-55[1] Gene H. Golub and Charles F. Van Loan, "Matrix Computations 4th Edition: Algorithm 9.3.1," JHU Press, 2013. |
Hi numo authors
I'm currently working with kernels and a couple of them needs the expm operation. in my project, I tried to implemented this operation using the Padé approximation using Nmatrix library. The implementation doesn't work (the error rate was huge) and several months after, I change to your library. Your library doesn't support the expm operation so I solve this problem using the Pycall library to use the expm operation from scipy.linalg. The problem is the huge consumption of memory and the time used copying the matrix ruby -> python -> ruby with very large matrix. Could you add the expm operation to your library?
I can give you as little reference my old code:
https://github.com/ElenaRojano/NetAnalyzer/blob/master/lib/NetAnalyzer/numo_expansion.rb
In line 79 is defined the expm operation, currently invoking python to solve it. But the commented code has a few implementations (with old Nmatrix code, but I thisk that you can guess easily the similar operations) that I tried to use with the cited resources. I think that the best reference is the scipy itself:
https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/matfuncs.py
The function is in the line 550.
Also, I would like to cite you in my papers. Do you have any citation that I could use?
Thank you very much in advance
Pedro Seoane
The text was updated successfully, but these errors were encountered: