Arls.py

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

(See http://www.rejtrix.net for similar software for other languages)

 

 

What this is:

Arls.py 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’ solvers and call arls:

from arls import arls1

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 arlsnn

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 arlseq

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 arlsgt

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

 

Case 5. You can combine all the above features:

from arls import arlsall

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

 

Case 6. When regularizing systems of equations, it is good to use as few and as unintrusive constraints as possible. For example, if constraining the first and last solution value to be >= 0.0 elicits a  good solution, then it is preferable to just do that, and nothing more. For another example, if you know that the solution has to be non-decreasing (or “rising”), or concave down, or some such then it is good to use that constraint first and not add other issues unless necessary. Here is a routine that can constrain the solution in 18 different combinations:

from arls import arlshape

nonneg = 0    # for example

slope = 1     # for example

curve = -1    # for example

x = arlshape(A, b, nonneg, slope, curve)

Here, set the three parameters like this:

      nonneg = 1 for a nonnegative solution.

      nonneg = 0 to skip this constraint.

      slope =  1 for a rising (or non-decreasing) solution.

      slope = -1 for a falling (or non-increasing) solution.

      slope =  0 to skip the slope constraint

      curve =  1 for a concave up (or possibly straight) solution

      curve = -1 for a concave down (or possibly straight) solution.

      curve =  0 to skip the concavity constraint.

 

For example, x = ArlsShape(A,b,0,1,-1) would demand a solution to  Ax=b which is rising and concave down.

 

Note: If all you want is nonnegativity, please use arlsnn(), which is more robust and faster.

 

Note: You may think it odd that we can add up to almost 3*m constraints to assure the shape is as you requested. That is way more constraints than variables! The reason this can make sense is that there is a natural linear dependence between these added equations such that only m of them can ever be invoked. What happens is that the other  approximately 2*m equations collapse to zeros and are then discarded by arls(). Further, note that it IS possible that a full “m” constraints MIGHT be invoked, leaving no place for the actual original equations. But typically this is the case only when the solution degenerates to all zeros. All zeros is always a possible answer because by our definition above, a nonnegative solution can be all zeros; and any flat line is considered an acceptable default case of either “rising” or “falling”, and a straight line (of any slope) is considered an acceptable default case of either “concave up” or “concave down”. So a flat line at 0.0 will always satisfy this routine’s constraints -- but maybe not YOUR desires!

 

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  rejones7@msn.com