A Python Module for Automatic Regularization of Ill-Conditioned Linear Systems

(See for similar software for other languages)



What this is: is a complete package for AUTOMATICALLY solving difficult linear algebraic systems (or “matrix problems”).  Basically, these methods completely free the user who has a difficult, usually ill-conditioned, problem from worrying about :

·         choosing a rank to use

·         choosing a level of singular value to zero when using a Generalized Inverse

·         estimating the error in the matrix coefficients

·         estimating the error in the right-hand side

·         choosing the Tikhonov regularization parameter, lambda

Arls is available on the Python Package Index, PyPI. So you should be able to install it in you Python installation with pip:

                 Pip install arls

How to Use arls. We assume you have some equations to solve, and can use Python. Put the coefficients in rows of matrix A, and put the right hand side values in a column, b. It does not matter whether the number of rows equals the number of unknowns, or is larger or smaller than the number of unknowns, as arls’s routines handle all cases.  Singularity or ill-conditioning are expected and considered normal behavior by these methods.

·         If the number of rows is less than or equal to the number of unknowns and the data is consistent and accurate, then the solution will typically be exact except for round-off error, as any good standard linear equation solver would do.

·         If the number of rows is larger than the number of unknowns, or the data is significantly in error, then a best fit, or “least-squares” solution will be found.

·         If the system of equations had serious problems such as unintended dependency between rows, or “singularities”, or ill-conditioning, then arls will usually be able to address those issues automatically and still return a reasonable solution. This is the interesting case for arls, and the reason for its existence. Arls was developed to automatically solve so-called “ill-conditioned” systems which often otherwise require special handling or are considered impossible to solve. Of course, when arls’ algorithms produce a solution of such a difficult problem the residual error (that is, the values in Ax-b, or the single value, norm(Ax-b)) will of necessity be larger than what the residual would have been if a standard least-squares solver had been able to handle the problem. So always check the result for reasonableness and appropriateness to your technical needs.

Case 1. Once you have created your matrix, A, and the corresponding right-hand-side, b, we have the problem

     A x = b

To solve it, import arls’s solvers and call arls:

from arls import arls1, arlseq, arlsgt, arlsnn, arlsall

x = arls1(A,b)[0]

Case 2. If your system of equations results in some negative values, and that is not acceptable, then just call our non-negative solver:

from arls import arls1, arlseq, arlsgt, arlsnn, arlsall

x = arlsnn(A,b)

Case 3. If you have a special situation, such as needing the solution to add to, say, 100.0 percent, or the first value to start at exactly 0., or any such exact requirement (not just more regular “least-squares” equations) then put those requirements as rows in a separate set of equations, Ex = f and then solve the dual system, by calling our “equality constraint” solver:

from arls import arls1, arlseq, arlsgt, arlsnn, arlsall

x = arlseq(A,b,E,f)

For example, if there are 5 unknowns, to require the solution to add to 100, use this:

                E = [ 1, 1, 1, 1, 1]

              f = [100]                       

Note: Depending on what you put in Ex = f, it may be impossible to satisfy all the equations. For example, you might accidentally say x[0] = 1., and also x[0] = 2. In such cases arlseq() will attempt to retain the equations that affect the result most.  You can check the status of the equality constraints by examining the residual:

                 R = E*x – f                 (where x is the solution returned by arlseq() )

Normally all elements of R will be zero. If R[3], for example, is non-zero, then equation 3 was incompatible with the other constraints.  You should try to resolve what is going on with these constraints.


Case 4. If you have some special need for certain unknowns to satisfy given bounds, then put those requirements in a separate set of equations, Gx >= h. Be sure to express them all in “greater than” form, such as x[3]>1.0. If you need to use a “less than” equation, such as x[0] < 2.0, then multiply the equation on both sides by -1 and you will get a “greater than” constraint instead of a “less than” constraint”:

-x[0] >= -2.0. 

Note that using this technique you can bound a value from below and above, such as

                X[3] >= 2

x[3] <= 4                 ( a.k.a.  -x[3] >=  -4 )               

If we assume there are, say, 5 unknowns, then these two constraints can be represented by

G = [  [0, 0, 0, 1, 0],                                (remember Python indexing starts at 0, not 1)

          [0, 0, 0, -1, 0]  ]                             (you have to negate “less than” equations)

h = [ 2,  -4 ]

Then you can solve this combined system by calling:

from arls import arls1, arlseq, arlsgt, arlsnn, arlaall

 x = arlsgt(A,b,G,h)


Case 5. Lastly, you can combine all the above features:

from arls import arls, arlseq, arlsgt, arlsnn, arlsall

x = arlsall(A,b,E,f,G,h)


Extra Information: In the case of arls1() you can get extra information

x, nr, ur, sigma, lambda = arls1(A,b)

In this case the extra parameters have these meanings when b contains only one column:

  nr : (int) is the numerical rank of the matrix A. This measure of rank ignores singular values that are zero or near zero.

  ur : (int) is the computed “usable rank”, which is a function of both A and b, and tells the solver at what point the error in A and b combine to ruin the use of more than ur rows of the orthogonalized version of the system of equations. In other words, ur is the breakpoint in the orthogonalized version of the system of equation between rows in which the error in b is acceptably small, and rows in which the error in b is uselessly large.

  sigma : (float) is the estimated estimated right-hand-side root-mean-square error. This means that our analysis reveals the actual error in b – a unique feature!

  lambda : (float) is the Tikhonov regularization parameter that goes with the estimated right hand side error, sigma.


Lastly: if you have a special need, such as a solver that incorporates two or more of the above four cases, or a need to minimize some measure other than the sum of the squares of the errors of the solution, or whatever, please send us a note. We may be able to provide a special solver.


Contact us:

Rondall E. Jones