After some testing of Saullo G. P. Castro's answer I've found
- a part that can be made more efficient (also can be found in Matthew Lewis' answer)
- a bug that causes the sign to be incorrect
I described both in the GIST of LuizFelippe mentioned in the comments of Saullo G. P. Castro's answer. However, the GIST seems to be inactive, so I decided to post an answer here given that I also do not have enough reputation for a comment yet.
Inefficient computation of the determinant of the factor L
There is a tiny improvement possible, which I found out during debugging and from the SuperLU documentation which states that the matrix L is unit lower triangular, i.e., its main diagonal is a vector of ones:
So in principle, it should be possible to drop all the terms involving L because the sign will always be +1.0 and the logarithm of the respective product will be 0.0.
Bug in sign computation
Since only the row permutations but not the column permutations were included in the proposed code, there was a 50% chance for a wrong sign in the determinant (since the permutation matrices have a determinant of either +1 or -1 which leaves a 50/50 chance for a product of two such matrices to be +1 and -1, respectively). Even though there are tests mentioned in the GIST, this failed for me the first time I ran the code.
The column permutations can be included in the exact same way as the row computations, so the fixed code is given by the following.
When the line marked with <-- š£ is uncommented, this yields the original approach where the column permutations are not considered.
### Imports ###
import numpy as np
from scipy.sparse import linalg as spla
### Functions ###
def sparse_slogdet_from_superlu(splu: spla.SuperLU) -> tuple[float, float]:
"""
Computes the sign and the logarithm of the determinant of a sparse matrix from its
SuperLU decomposition.
References
----------
This function is based on the following GIST and its discussion:
https://gist.github.com/luizfelippesr/5965a536d202b913beda9878a2f8ef3e
"""
### Auxiliary Function ###
def minimumSwaps(arr: np.ndarray):
"""
Minimum number of swaps needed to order a permutation array.
"""
# from https://www.thepoorcoder.com/hackerrank-minimum-swaps-2-solution/
a = dict(enumerate(arr))
b = {v: k for k, v in a.items()}
count = 0
for i in a:
x = a[i]
if x != i:
y = b[i]
a[y] = x
b[x] = y
count += 1
return count
### Main Part ###
# the logarithm of the determinant is the sum of the logarithms of the diagonal
# elements of the LU decomposition, but since L is unit lower triangular, only the
# diagonal elements of U are considered
diagU = splu.U.diagonal()
logabsdet = np.log(np.abs(diagU)).sum()
# then, the sign is determined from the diagonal elements of U as well as the row
# and column permutations
# NOTE: odd number of negative elements/swaps leads to a negative sign
fact_sign = -1 if np.count_nonzero(diagU < 0.0) % 2 == 1 else 1
row_sign = -1 if minimumSwaps(splu.perm_r) % 2 == 1 else 1
col_sign = -1 if minimumSwaps(splu.perm_c) % 2 == 1 else 1
# col_sign = 1 # <-- š£ If this is uncommented, this produces the `perm_r`-only code
sign = -1.0 if fact_sign * row_sign * col_sign < 0 else 1.0
return sign, logabsdet
I implemented a more extensive test against numpy.linalg.slogdet (takes 5 to 10 minutes on an M4 MacBook Pro).
It tests at least 10 matrices for every given row/column count between 50 and 1000 to ensure consistency and not just lucky shots. Since we do not want to test SuperLU's ability to solve random sparse matrices which can be ill-conditioned, a matrix that cannot be solved will be regenerated in a random fashion.
While this test passes ā
with the suggested fix (the line with <-- š£ is left commented), it fails ā on the first attempt when using the original code (the line with <-- š£ is active).
### Tests ###
if __name__ == "__main__":
# Imports
import numpy as np
import scipy.sparse as sprs
from scipy.sparse.linalg import splu as splu_factorize
from tqdm import tqdm
# Setup of a test with random matrices
np.random.seed(42)
# n_rows = np.random.randint(low=10, high=1_001, size=20)
density = 0.5 # chosen to have a high probability of a solvable system
n_rows = np.arange(50, 1001, dtype=np.int64)
# Running the tests in a loop
for index in tqdm(range(0, n_rows.size)):
m = n_rows[index]
num_tests_passed = 0
num_attempts = 0
failed = False
while num_tests_passed < 10:
# a random matrix is generated and if the LU decomposition fails, the
# test is repeated (this test is not there to test the LU decomposition)
num_attempts += 1
matrix = sprs.random(m=m, n=m, density=density, format="csc")
try:
splu = splu_factorize(matrix)
except RuntimeError:
tqdm.write(
f"Could not factorize matrix with shape {m}x{m} and density "
f"{density}"
)
if num_attempts >= 100:
tqdm.write(
f"Could not generate a solvable system for matrix with shape "
f"{m}x{m}"
)
failed = True
break
continue
# first, the utility function is used to compute the sign and the log
# determinant of the matrix
sign, logabsdet = sparse_slogdet_from_superlu(splu=splu)
# then, the sign and the log determinant are computed by NumPy's dense
# log determinant function for comparison
sign_ref, logabsdet_ref = np.linalg.slogdet(matrix.toarray())
# the results are compared and if they differ, the test is stopped
# with a diagnostic message
if not (
np.isclose(sign, sign_ref) and np.isclose(logabsdet, logabsdet_ref)
):
print(
f"Failed for matrix with shape {m}x{m}: "
f"sign: {sign} vs. {sign_ref} and "
f"logabsdet: {logabsdet} vs. {logabsdet_ref}"
)
failed = True
break
# if the test is successful, the loop is continued with the next iteration
del splu
num_tests_passed += 1
if failed:
break
