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 1).
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.
1) 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.
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
pentapysolves the systemAx = b, it basically factorisesAand transformsbas follows for algorithm I:A = LUwhereLis a lower triangular matrix whose entries will not be discussed here (its lower triangle is dense) andUis the unit upper triangular matrix whereSince the main diagonal of
Uholds 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 ofAis completely encoded inLfor which we then also need the main diagonal elements given that this triangular as well. Consequently, the same computations as forUapply. 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) @ bWhat does that help us? Well, we now need to look at how the update withinv(L)is achieved:The parts I marked in red help us getting the main diagonal of
inv(L)because this means thatCombining 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 1).
I tested this for algorithm I against
np.linalg.slogdetand it works ✅.If we allow the user to set a flag like
with_logdet, this could be returned along or not. If the default isFalse, 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 usepentapyto a way stronger extent than we currently do.1) 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.