neo

neo copied to clipboard

The computation of eigenvalues can be used for the determinant, or perhaps there is a more specific LAPACK function

The Faq mentions:

Are there routines in LAPACK to compute determinants?

No. There are no routines in LAPACK to compute determinants. This is discussed in the "Accuracy and Stability" chapter in the LAPACK Users' Guide.

In gonum this is done using the LU factorization of the matrix.

var lu LU
lu.Factorize(a)
lu.LogDet()
det, sign := lu.LogDet()
res := math.Exp(det) * sign

Yeah, I know it's not very stable. There is no LU factorization here yet, but the Schur form can be used for now (although it is more costly; the determinant will be unstable for large matrices anyway)

https://github.com/unicredit/neo/commit/742c60f852a4d3ebc531b7089fb4cb9d2068b916 introduced the two functions. det can be optimized, but as remarked it only makes sense for small matrices anyway, so it not much of an issue for now

I'm too new to nim and this to make a PR, but I cobbled together a cholesky decomp. I'm not sure if I used the T right since I think that has to be float64 with the dpotrf call.. anyway here it is:

proc cholesky*[T](a: Matrix[T]): Matrix[T] =
  var
    h = a.clone()
    n = a.ld.cint
    info: cint
    ulo: cstring
  ulo = "L"
  dpotrf(ulo, addr n, h.fp, addr n, addr info)
  if info > 0:
    raise newException(LinearAlgebraError, "ERROR finding the LU decomp")
  for i in 0 ..< n:
    for j in i + 1 ..< n:
      h[i, j] = 0.0
  return h

awesome work on this package!