Numerical solving system of equations (ill-conditioned)

How to solve a system of equations Ax=b when the coefficient matrix is ill-conditioned? Such a matrix is almost singular, and the computation of its inverse, or solution of a linear system of equations is prone to large numerical errors.

There are several ways to estimate or calculate the condition number of a matrix (here)but just consider that condition number of a singular matrix is infinity and the best condition number for solving Ax=b is 1.

What happen if you try to solve such a system of equations?
We need a matrix with high condition number, consider the Hilbert matrix defined by: A [N by N] (aij = 1/(i+j-1)). The right hand side vector (b) will be considered as b = (1, 1, … , 1).
We use Hilbert matrices of various orders as test problems.

Hilber Matrix

If you try to solve Hilber matrix with different orders you easily loose precision by increasing the condition number.

Rule of thumb

If the linear system is solved in m-digit floating point arithmetic
and if the condition number of A is of the order 10^α , then only
m − α − 1 digits in the solution can be trusted.

What if we have matrix with condition number of 10²⁸ or higher?
As you can guess there will be no accurate results using standard precision numbers with 16 decimal digits. Numpy also provide long double (numpy.float128 with 33 decimal digits) but the problem is that (As far as I know) none of the standard solver methods in Scipy and Numpy support float128.

Solution

We need to use multiprecition libraries and do the calculation in higher precision. For example if we need at least 10 meaningful decimal digits with condition number = 10³⁰, we need to consider 41 decimal digits.

The MPFR library is a C library for multiple-precision floating-point computations with correct rounding.

Eigen support MPFR multi-precision data structure, so we can simply define matrix with arbitrary number of decimal digits and use default solvers for dense and sparse matrices.

using MatrixXmp = Matrix<mpreal, Dynamic, Dynamic>;
using VectorXmp = Matrix<mpreal, Dynamic, 1>;
mpfr::mpreal::set_default_prec(mpfr::digits2bits(num_decimal));
mpfr::digits2bits function calculate the required number of bits for number of decimals.
for (int n = 4; n < nn; ++n)
{
 MatrixXmp A = hilbert_matrix(n);
 VectorXmp x(n);
 for (int i = 0; i < n; ++i)
 x(i) = 1.;
 VectorXmp b = A * x;
 VectorXmp xx = A.partialPivLu().solve(b);
}

Using Eigen and MPFR to solve ill-conditioned Ax=b.

The full example for this post is available here.
To compile the code use g++ example.cpp -lmpfr -lgmp and make sure to install MPFR and GMP libraries.