Julia not solving sparse Linear system

alexp9 0
Viewed 75

Introduction

I'm doing research in computational contact mechanics, in which I try to solve a PDE using a finite difference method. Long story short, I need to solve a linear system like Ax = b.

The suspects

In the problem, the matrix A is sparse, and so I defined it accordingly. On the other hand, both x and b are dense arrays.

In fact, x is defined as x = A\b, the potential solution of the problem.

So, the least one might expect from this solution is to satisfy that Ax is close to b in some sense. Great is my surprise when I find that

julia> norm(A*x-b) # Frobenius or 2-norm
5018.901093242197

The vector x does not solve the system! I've tried a lot of tricks discover what is going on, but no clues as of now. My first candidate is that I've found a bug, however I need more evidence to make this assertion.

The hints

Here are some tests that I've done to try to pinpoint the error

  • If you convert A to dense, the solution changes completely, and in fact it returns the correct solution.
  • I have repeated the proccess above in matlab, and it seems to work well with both sparse and dense matrices (that is, the sparse version does not agree with that of Julia's)
  • Not all sparse matrices cause a problem. I have tried other initial conditions and the solver seems to work quite well. I am not able to predict what property of the matrix can be causing this discrepancy. However;
  • A has a condition number of 120848.06, which is quite high, although matlab doesn't seem to complain. Also, the absolute error of the solution to the real solution is huge.

How to reproduce this "bug"

  1. Download the .csv files in the following link
  2. Run the following code in the folder of the files (install the packages if necessary
using DelimitedFiles, LinearAlgebra, SparseArrays;
A = readdlm("A.csv", ',');
b = readdlm("b.csv", ',');
x = readdlm("x.csv", ',');

A_sparse = sparse(A);

println(norm(A_sparse\b - x)); # You should get something close to zero, x is the solution of the sparse implementation

println(norm(A_sparse*x - b)); # You should get something not close to zero, something is not working!

Final words

It might easily be the case that I'm missing something. Are there any other implementations apart from the usual A\b to test against?

julia
1 Answers

To solve a sparse square system Julia chooses to do a sparse LU decomposition. For the specific matrix A in the question, this decomposition is numerically ill-conditioned. This is evidenced by the cond(lu(A_sparse).U) == 2.879548971708896e64. This causes the solve routine to make numerical errors in turn.

A quick solution is to use a QR decomposition instead, by running x = qr(A_sparse)\b.

The solve or LU routines might need to be fixed to handle this case, or at least maintainers need to know of this issue, so opening an issue on the github repo might be good.

(this is a rewrite of my comment on question)

Dan Getz0
Related