ruby-numo / numo-linalg

Linear Algebra Library for Ruby/Numo::NArray
BSD 3-Clause "New" or "Revised" License
38 stars 9 forks source link

Fix errors in the eigh method and to enable to solve generalized eigenvalue problem #9

Closed yoshoku closed 6 years ago

yoshoku commented 6 years ago

Problem

Numo::Linalg.eigh method returns the array consists of eigenvectors and eigenvalues. In the document, the first element of the returned array is the eigenvalues and the second element is the eigenvectors, but in fact the first element is eigenvectors. Moreover, when I use the eigh method with setting turbo option to true, Numo::NArray::CastError occurs;

> a = Numo::DFloat.new(3, 3).rand
> a = 0.5 * (a + a.transpose)

> Numo::Linalg.eigh(a)
=> [Numo::DFloat#shape=[3,3]
[[-0.572452, -0.70266, 0.422573],
 [0.639884, -0.0605972, 0.766078],
 [-0.512686, 0.70894, 0.48431]], Numo::DFloat#shape=[3]
[-0.410378, -0.0231387, 1.10204]]

> Numo::Linalg.eigh(a, turbo: true)
Numo::NArray::CastError: cannot cast to Numo::DFloat
...

Solution

I found the cause of these errors is that syev and sygv (hegv and heev) functions are used together in the eigh method.

case blas_char(a)
when /c|z/
  func = turbo ? :hegv : :heev
else
  func = turbo ? :sygv : :syev
end

Thus, I fixed to use only the sygv (hegv) and sygvd (hegvd) functions in the eigh method. The sygvd and hegvd functions are used for the turbo option. In addition, since these functions solve the general eigenvalue problem, I made the eigh method allow to specify the second matrix optionally like scipy.linalg.eigh.

Test

By running the following code, I checked the changed codes with this pull request.

https://gist.github.com/yoshoku/b2ecf1698a59c69f00fc79a34bc5926e