ruby-numo / numo-linalg

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

[Propose] Add new method for Cholesky factorization #19

Closed yoshoku closed 6 years ago

yoshoku commented 6 years ago

I would like to add new method for Cholesky factorization like scipy.linalg.cholesky.

Cholesky factorization decomposes a symmetric/Hermitian positive define matrix into the product of an upper / lower triangular matrix such as A = U^T * U. The already existing cho_fact method returns the matrix which has the Cholesky factor in upper / lower triangular part. However, the remain part consists of random values. Thus, it is difficult to reconstruct the original matrix using the returned matrix from the cho_fact method.

> a = Numo::DFloat.new(3,3).rand
> a = a.dot(a.transpose)
=> Numo::DFloat#shape=[3,3]
[[1.22905, 0.817261, 0.924107],
 [0.817261, 0.881999, 0.448001],
 [0.924107, 0.448001, 0.830497]]

> u = Numo::Linalg.cho_fact(a)
=> Numo::DFloat#shape=[3,3]
[[1.10862, 0.737185, 0.833562],
 [0.817261, 0.581856, -0.286134],
 [0.924107, 0.448001, 0.231946]]
irb(main):031:0> u.transpose.dot(u)
=> Numo::DFloat#shape=[3,3]
[[2.75094, 1.70679, 0.904603],
 [1.70679, 1.0827, 0.551913],
 [0.904603, 0.551913, 0.830497]]
irb(main):032:0>

In contrast, the cholesky method in this PR returns the upper / lower triangular matrix consits of Cholesky factor. The original matrix can be reconstructed by the product of the matrix returned by the cholesky method.

> a
=> Numo::DFloat#shape=[3,3]
[[1.22905, 0.817261, 0.924107],
 [0.817261, 0.881999, 0.448001],
 [0.924107, 0.448001, 0.830497]]

> u = Numo::Linalg.cholesky(a)
=> Numo::DFloat#shape=[3,3]
[[1.10862, 0.737185, 0.833562],
 [0, 0.581856, -0.286134],
 [0, 0, 0.231946]]
> u.transpose.dot(u)
=> Numo::DFloat#shape=[3,3]
[[1.22905, 0.817261, 0.924107],
 [0.817261, 0.881999, 0.448001],
 [0.924107, 0.448001, 0.830497]]