79408928

Date: 2025-02-03 13:17:58
Score: 1
Natty:
Report link

After some testing of Saullo G. P. Castro's answer I've found

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: image

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
Reasons:
  • RegEx Blacklisted phrase (1.5): do not have enough reputation
  • Long answer (-1):
  • Has code block (-0.5):
  • Low reputation (1):
Posted by: MothNik