[Enhancement] Evaluate log-determinant for solve

[MothNik]

While #27 is open and changes a lot of internal logic, I want to propose another feature for a version upgrade: the computation of the log-determinant. It would add into #27 very easily with a few lines of code.

If we look at how pentapy solves the system Ax = b, it basically factorises A and transforms b as follows for algorithm I:

• A = LU where L is a lower triangular matrix whose entries will not be discussed here (its lower triangle is dense) and U is the unit upper triangular matrix where

  • the main diagonal consists of ones,
  • the first superdiagonal consists of the $\alpha_{i}$-values, and
  • the second superdiagonal consists of the $\beta_{i}$-values
  • all the rest is zeros

  Since the main diagonal of U holds only ones and the log determinant of a triangular matrix is simply the sum of the logarithms of its main diagonal elements, this is 0.
  So, the determinant of A is completely encoded in L for which we then also need the main diagonal elements given that this triangular as well. Consequently, the same computations as for U apply. But how should this be done? It cannot be computed when it's dense because that would ruin performance. Well, this is where the next point comes in:

• zeta = inv(L) @ b
  What does that help us? Well, we now need to look at how the update with inv(L) is achieved:

  

  The parts I marked in red help us getting the main diagonal of inv(L) because this means that

  

Combining all of this with the fact that the main diagonal elements of the inverse of a triangular matrix are just the reciprocals of its main diagonal elements, we get

Long story short, that means that if we sum up the logarithms of the $\mu_{i}$-values that we get from the factorization anyway (keeping care of the sign), we can get the log-determinant in no time at all. The same applies for $\psi_{i}$ for algorithm two, so we get:

For multiple right-hand sides, the additional load will be vanishing because the factorization and determinant computation only need to be done once.
The error case of $ln\left(0\right)$ is already caught for both algorithms I and II because the factorization will throw a zero-division error before any logarithm is applied ¹⁾.

I tested this for algorithm I against np.linalg.slogdet and it works ✅.

If we allow the user to set a flag like with_logdet, this could be returned along or not. If the default is False, the Python -API does not change, but we would need to ⚠️ adapt the C-level-API again by adding an additional bool and a double pointer ⚠️

In chemotools, we rely heavily on the log determinants for Bayesian regression, but currently we use SciPy's solvers for that. Therefore, this feature would enable us to use pentapy to a way stronger extent than we currently do.

¹⁾ The determinant can however NOT be used to check whether the matrix is singular. This is a common mistake that comes from math classes where everything was done in 100% exact arithmetics. On a computer with finite precision, this is the worst check to be performed and only the condition number can provide true insights for this.