Arls.go

A “Go” Package for Automatic Regularization of

Ill-Conditioned Linear Systems

 

(See http://www.rejtrix.net for similar software for Python and C++.)

 

What this is:

Arls.go is a complete package for AUTOMATICALLY solving difficult linear algebraic systems (or “matrix problems”).  These algorithms have been available in C++ for several years, but we recently decided to migrate them to Python and Go. Please see the web pages for those two packages at http://www.rejones7.net/ for examples, tutorials, etc. , as they mostly apply for Go, except for their extra, more specialized solvers that accept inequality and equality constraints.

Basically, all these three packages completely free the user who has a difficult, usually ill-conditioned, matrix 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

They work equally well for well-behaved systems too, so you can use them for all your dense matrix problems. (The execution time will be longer than most other solvers.)

You may download our Go package here:  arls.go.

How to Use Arls. We assume you have some equations to solve, Ax = b.

 Create A as a *mat.Dense, and b as a  *mat.VecDense. Put the coefficients in rows of matrix A, and put the right hand side values in b. It does not matter whether the number of rows, m, equals the number of unknowns, n, 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.

Case 1. To solve Ax =b for x, import Arls’s solvers and Go’s “math” and “mat” packages, and call Arls:

import ( "math"; "gonum.org/v1/gonum/mat"; "arls")

A = …whatever…

b = …whatever…

x, nr, ur, sigma, lambda := Arls(A,b)

The return values are:

               x :(*mat.VecDense) is the solution.

nr :(int) is the traditional numerical rank of the matrix A. This measure of rank excludes singular values that are 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 :(float64) is the estimated right-hand-side root-mean-square error. This means that our analysis computes a useful estimate of thevactual error in b!

lambda :(float64) is the resulting Tikhonov regularization parameter.

Some details:

·         If the number of rows is less than or equal to the number of unknowns, and if the data is consistent and accurate then the solution will typically be almost exactly correct.

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

·         If the system of equations has 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.

·         Generally, for best results with Arls you should scale the rows of A such that the right-hand-side values (in b) have about the same estimated error. If you have no knowledge of the error in b, then probably you should scale the rows so that the sizes of the coefficients in A do not vary greatly between rows.

Case 2. If you need to solve many systems of equations use the SAME MATRIX, A, then you can save a great deal of execution time by computing the Singular Value Decomposition of A, which Arls must have, and reusing it. Here the calling sequence is Arlsvd(svd mat.SVD, b *mat.VecDense), and the returned parameters as exactly as for Arls().

import ("math"; "gonum.org/v1/gonum/mat"; "arls")

A = …whatever…

var svd mat.SVD

ok := svd.Factorize(A, mat.SVDThin)

if !ok { panic("SVD failed to factorize A") }

b1 = …whatever…

x1, nr1, ur1, sigma1, lambda1 := Arlsvd(svd, b1)

b2 = …whatever…

x2, nr2, ur2, sigma2, lambda2 := Arlsvd(svd, b2)

      . . .

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

import ( "math"; "gonum.org/v1/gonum/mat"; "arls")

A = …whatever…

B = …whatever…

x, nr, ur, sig, lam := Arlsnn(A,b)

There is no version of Arlsnn that accepts the SVD of A, as that is not really helpful in this case.

Case 4. Sometimes it is important that certain equations be given priority over others. In the case of systems being solved by least-squares, you can somethings scale that row of the system to much larger values than the other rows as a way to give that row (or rows) priority. But it is also possible to mathematically adjust the system so that some equations are solved exactly, and the others solved as usual. (Of course, there is no magic here… if your “exact” equations (a.k.a. “equality constraints”) conflict with each other, Arls will of necessity eliminate certain of the equations until a consistent set is found. This process is complex and not fool-proof.) So you will need to split the equations in to two set:

                A x = b                           for the regular (least-squares) equations

                E x ==  f                                 for the exact, or “equality” equations.

Then call

x, nr, ur, sig, lam := Arlseq(A,b,E,f)

to solve these equations. The return values are as for Arls.

Case 5. Sometimes you may know certain limits on the solution values. For example, above we discussed the need to force all solution values to be non-negative. But you may also know that the solution should never exceed 1. To use our “inequality constraint” solver you must form all these conditions into equations with a “greater than” symbol. (If you need “less than” then you must negate both sides of the equation to convert it to a “greater than” equation.)

For example, one might have

                Xi <= 1.0  for all i

 or           X1 >= 1, x2>=2, and x9>=5

or            X1 + x2 + x3 + … +x10 >= 100.

Then put these equations into matrix form, like this:

     Gx >= h

and call

x, nr, ur, sig, lam := Arlsgt(A,b,G,h)

Case 6. Finally, you are allowed to use all these features at once. That is if you have

                A x = b                    for the regular (least-squares) equations

                E x == f                   for the exact, or “equality” equations.   

G x >= h                        for the “greater than” equation

Then you can call

x, nr, ur, sig, lam := Arlsall(A,b,E,f,G,h).

 

Contact us:

Rondall E. Jones,

rejones7@msn.com

Sandia Laboratories, Albuquerque, 1967-2011

Publications: Please Google “Rondall Jones Sandia” and ignore results not containing the exact spelling of “Rondall”

Ph.D., University of New Mexico, 1985.  Dissertation: Solving Linear Algebraic Systems Arising in the Solution of Integral Equations of the First Kind; Advisor: Cleve B. Moler, creator of MATLAB and co-founder of MathWorks. My thanks to Richard Hanson (deceased), co-author of Solving Least Squares Problems, for additional guidance for this work.