UMFPACK: enforcing unsymmetric strategy with symmetric matrix may yield wrong results

[cpmech]

Describe the bug

Note: this is an experiment with UMFPACK_STRATEGY_UNSYMMETRIC, and may not be a bug, really.

Tested on Ubuntu 22.04 and Ubuntu 23.04 with UMFPACK V5.7 (Oct 20, 2019) and UMFPACK V6.2 (Sept 18, 2023).

When running UMFPACK on two similar symmetric matrices, the error may increase significantly if we enforce the unsymmetric strategy (clearly not recommended as stated in the documentation).

Steps:

1. Generate two sparse and symmetric matrices for a simple (linear elastic) finite element analysis. One matrix (mat # 2) corresponds to a slightly finer mesh than the other (mat # 1).
2. Implement UMFPACK to solve a linear system with those matrices and a unit right-hand side.
3. Check the behavior of UMFPACK with UMFPACK_STRATEGY_UNSYMMETRIC and without it (Automatic); i.e., enforce the unsymmetric strategy or not.
4. Observe the following:
5. With the first matrix (mat # 1), by enforcing the unsymmetric strategy, the residuals (error) increase from 10^-12 to 10^-5
6. With the second matrix (mat # 2), by enforcing the unsymmetric strategy, the residuals (error) increase from 10^-12 to 10^+3 (too high).

To Reproduce

Use the modified version of umfpack_di_demo.c that comes with SuiteSparse. The modified version is available here: https://github.com/cpmech/umfpack-matrix-market (look for solve_matrix_market.cpp). This modified version first loads a matrix written as MatrixMarket.

After installing libsuitesparse-dev on Ubuntu, the results can be obtained by running bash ./compare-1-and-2.bash.

Alternatively, the problem may be directly reproduced using Docker:

docker pull cpmech/umfpack-matrix-market
docker run --rm -it cpmech/umfpack-matrix-market:latest /bin/bash
bash ./compare-1-and-2.bash

Expected behavior

If this is a bug, UMFPACK should report that the matrix is best solved as symmetric.

If this is not a but, the documentation (which already is excellent) should strongly (even more strongly) discourage setting UMFPACK_STRATEGY_UNSYMMETRIC.

Screenshots

... FILE: pressurized-cylinder-linear-elastic-symmetric-1.mtx
... UMFPACK V6.2 (Sept 18, 2023) demo: _di_ version
... SUCCESS: matrix loaded
... SUCCESS: CSC arrays allocated
... SUCCESS: COO converted to CSC
... SUCCESS: symbolic factorization completed
... SUCCESS: numeric factorization completed
... SUCCESS: solution calculated
... max_norm of residual: 5.00222e-12
... max_diff using COO  : 5.00222e-12
... SUCCESS: numerical solution is within tolerance

... FILE: pressurized-cylinder-linear-elastic-symmetric-1.mtx
... ENFORCING UNSYMMETRIC STRATEGY
... UMFPACK V6.2 (Sept 18, 2023) demo: _di_ version
... SUCCESS: matrix loaded
... SUCCESS: CSC arrays allocated
... SUCCESS: COO converted to CSC
... SUCCESS: symbolic factorization completed
... SUCCESS: numeric factorization completed
... SUCCESS: solution calculated
... max_norm of residual: 1.09071e-05
... max_diff using COO  : 1.09071e-05
... FAIL: the error is too high

... FILE: pressurized-cylinder-linear-elastic-symmetric-2.mtx
... UMFPACK V6.2 (Sept 18, 2023) demo: _di_ version
... SUCCESS: matrix loaded
... SUCCESS: CSC arrays allocated
... SUCCESS: COO converted to CSC
... SUCCESS: symbolic factorization completed
... SUCCESS: numeric factorization completed
... SUCCESS: solution calculated
... max_norm of residual: 5.05906e-12
... max_diff using COO  : 5.05906e-12
... SUCCESS: numerical solution is within tolerance

... FILE: pressurized-cylinder-linear-elastic-symmetric-2.mtx
... ENFORCING UNSYMMETRIC STRATEGY
... UMFPACK V6.2 (Sept 18, 2023) demo: _di_ version
... SUCCESS: matrix loaded
... SUCCESS: CSC arrays allocated
... SUCCESS: COO converted to CSC
... SUCCESS: symbolic factorization completed
... SUCCESS: numeric factorization completed
... SUCCESS: solution calculated
... max_norm of residual: 4189.11
... max_diff using COO  : 4189.11
... FAIL: the error is too high

Desktop (please complete the following information):

For the local test code:

• OS: Ubuntu 23.04 Linux
• compiler: g++ (Ubuntu 12.3.0-1ubuntu1~23.04) 12.3.0

From the Docker image/container:

-OS: Ubuntu 22.04.3 LTS -compiler: g++ (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0

Additional context

Nope. Thanks and best regards.

[DrTimothyAldenDavis]

I'll give it a try with your particular matrix once I get back to the office (recovering from surgery at the moment), but a few thoughts until then.

This isn't likely to be a bug. I've seen it happen before in other matrices and it's a result of an interplay between the NP hard problem of trying to find a good sparsity preserving ordering and one that is also numerically acceptable. It's impossible to find a good heuristic that satisfies both. I've seen some symmetric matrices where unsymmetric pivoting can lead to catastrophic fillin and numerical errors. But there are also many cases where the unsymmetric strategy is vastly superior in terms of fillin, memory etc.

So there's no way to provide a firm answer as to what strategies will work and which may fail. I can't say "never use the unsymmetric strategy", because my heuristics for deciding which strategy to use cannot be perfect. I might guess wrong. The symbolic ordering is an NP hard problem after all. Toss in floating point round off issues and it becomes even harder.

QR factorization is more reliable in terms of accuracy so it doesn't have this numerical problem. However, the memory and time taken can be so large as to fail on those grounds, by running out of memory and time. It too has to use an ordering heuristic for an NP hard problem.

So it would not be reasonable to revise the documentation by much. I will check (later...) to see if I mention the possibile impact of the ordering strategy on numerical accuracy in the docs, and add a note to that effect. I might already say that; if not I should do so.

[cpmech]

Good morning! Thanks for looking into this.

Below are some screenshots comparing the UMFPACK infos for the 1st ("good") and the 2nd ("bad") matrices.

With AUTO mode enabled

Left = "Good" matrix; Right = Bad matrix

With UNSYMMETRIC mode enabled

Left = "Good" matrix; Right = Bad matrix

Observations

The AUTO mode results in the following rcond (estimated reciprocal condition number):

1. Good: rcond = 1.54e-02
2. Bad: rcond = 2.60e-02

The UNSYMMETRIC mode results in the following rcond (estimated reciprocal condition number):

1. Good: rcond = 4.31e-13
2. Bad: rcond = 5.91e-17 (aka, "singular", that's probably why the solution A*x=b fails)

I'm still wondering why the "simple" FEM stiffness matrix is causing these results. The mechanical problem I'm solving is a 2D linear elastic plane-strain ring representing a pressurized cylinder (I'll generate a figure of the mesh later).

Cheers. Dorival

[DrTimothyAldenDavis]

I've added a comment about the umfpack strategy, for the v6.2.1 release of UMFPACK.

[cpmech]

Thank you and best wishes ;-)