[SciPy-Dev] Resubmission of autosol.py

[rondall jones]

I am resubmitting an enhanced version of my autosol.py module. I have uploaded the code to scipy\linalg, and have uploaded a silent full-coverage test to scipy\linalg\tests.
WHAT IS IN THIS MODULE
Autosol.py is a package of automatically regularized solvers for dense linear algebraic systems. The term “regularize” (see Wikipedia) means to add extra information to a set of equations (or other application) in order to move a poor result more toward what is needed. In our case this extra information is mainly three things: 1) requesting the solution to exclude abnormal behavior caused by unavoidable noise in the data; 2) adding to the system of regular least-squares equations some specific equations that are known to be exactly true, such as “the solution sums to 1.0”; 3) giving the solver specific bounds on the answer by way of inequality constraints, such as requesting a non-negative solution.
arls() is the main automatically regularized least squares solver for Ax =b. The equations can be easy to solve, mildly or strongly ill-conditioned, or singular. “arls” stands for “Automatically Regularized Least Squares.” I have changed the names of the routines in autosol from the original long names to more scipy-like short names. I have left the module name as it was.
arlseq() builds on arls() to add Equality Constraints. There can be any number of equality constraints. Conflicts or other issues between these equations are resolved. Each equation is either solved exactly or rejected as incompatible with the other equations. This is normal behavior for such solvers.
arlsgt() builds on arlseq() to add Inequality Constraints. These often have the form of “the sum of the solution elements must be 100”, or “the solution must start at 1.0”, etc. There can be any number of inequality constraints. Inequality constraints are promoted to equality constraints when they are violated. This is, again, normal behavior for such solvers. See Lawson & Hansons’ “Solving Least Squares Problems” book.
arlsnn() builds on arlsgt() to provide a simple non-negativity constraint, like scipy’s nnls().
WHAT IS SIGNIFICANT ABOUT THESE ROUTINES?
These solvers are different than anything currently in scipy. (Skip this paragraph if you are not interested in the mathematical detais!) Briefly, if Ax = b, and A = USVt (an SVD) then
USVt * x = b can be written S(Vt*x)= Ut*b. The “Discrete Picard Condition” that is the basic insight that autosol uses says that the right-hand side, Ut*b, must decrease faster than the singular values in S. In ill-conditioned problems that condition fails, often dramatically, and of course it fails by definition if A is singular. In these cases a careful analysis of the vector g = (S+ * Ut*b) can produce what might be call a “usable rank” (which is smaller than the numerical rank) which then allows us to produce an estimate of the error in b, which then directly leads to an appropriate value for the lambda regularization parameter used in Tikhonov regularization. Thus, we can regularize Ax=b with no user input hints, and with no failure-prone iterative calculations. This method is very robust. It can handle multiple difficulties such as linear dependencies between rows, small singular values, zero singular values, high error levels in b, etc., etc., with no danger whatsoever that a numeric overflow or other error will occur. If the SVD does not crash (and such crashes are very rare indeed) then autosol’s algorithms will complete normally. There are no error modes to report. And there is no need for any special user input to guide the process or limit iteration counts, etc.
The other solvers – arlseq(), arlsgt(), and arlsnn() build in a somewhat classic fashion on arls(), using traditional methods, though those methods have been generalized or enhanced in places, such as how problematic Equality Constraints are rejected.
HOW DOES THIS PACKAGE HELP SCIPY?
In contrast to autosol, all current solvers in scipy work “blindly”. That is, they work without knowing whether the Ax=b problem they are solving is actually ill-conditioned or not. For example, lstsq() only diagnoses a problem when A is fully singular. But even a moderately ill-conditioned system will cause lstsq() to deliver a bad, high-norm, often oscillatory solution without notifying the user. This happens because the singular values only have to be small compared to the error level in b – not near zero -- to produce a pathological result.
On the other hand, lsmr(), which is the current primary routine for handling ill-conditioned systems (and often produces beautiful results) is also “blind” in that it never knows for sure whether the problem to which it is applying a regularization algorithm is actually ill-conditioned. The result is that it will often provide an inappropriately smoothed answer to a perfectly easy problem. Also, lsmr() is prone to producing erratic results with minor changes to a problem. I have prepared a brief illustration of these problems here: http://www.rejones7.net/Autosol/DemoArlsVsLsmr/index.htm
Additionally, the non-negative solver provided in autosol has advantages over the classis nnls() in scipy. The classic nnls() algorithm (from Lawson & Hanson’s classic Solving Least Squares Problems) works very hard at producing the smallest possible residual, and often ends with a failure to converge. It also causes problems with zeroing an excessive number of columns. For example, in a recent consultation regarding a biochemical problem, nnls() on average deleted over 50% of the columns in the problem to produce a non-negative result. But arlsnn() on average deleted of only 17% to achieve non-negativity, thus remaining more relevant to the model. And the slightly higher residual obtained by arlsnn() is usually trivial. Also, arlsnn() cannot fail to converge.
WHAT ISSUES ARE THERE?
One question involves computational cost. I have timed a number of problems being solved with lsmr() and autosol’s arls(). The timing is almost identical. On the other hand, lstsq() is greatly faster than lsmr() and arls(), so there is a cost for regularization, whether with lsmr() or arls(). The advantage of autosol’s algorithm is that it knows what difficulties it is dealing with: it handles well-conditioned problems with full accuracy, like lstsq(), and handles ill-conditioned problems as well as lsmr(), without ever foolishly applying a regularization technique to a well-conditioned problem.
Of course, no algorithm or software package is without flaws. Like lsmr() -- but less often than lsmr() -- autosol’s algorithms will fail to see ill-conditioning where some exists, or “see” ill-conditioning that is not really present. We have tuned the detection method in the primary solver, arls(), for a decade, in a C++ environment, and have enhanced its reliability significantly during redevelopment of these routines for Python. It should be significantly more reliable then lsmr(), according to the behavior I have seen in testing.
Some history: I developed the first version of the key algorithms in this package as part of a Ph.D. dissertation in the 1980’s. Due to a career change I did not pursue further development of the methods far a long time. But I eventually improved the algorithms and moved the code from Fortran to C++ and made the resulting codes available on my web site for a decade or so. During this time the web site drew interesting customers from around the world, whose problem proved the usefulness of the heuristics in the arls() algorithm. Recently I decided to migrate the code to Python -- which has been fascinating -- and resulted in some further refinements. See http://www.rejones7.net/ for both the C++ material and this Python material and related information and examples.
I am hoping that the scipy community will accept this package – with any necessary changes that are required by scipy – and let users learn about it as an alternative to the more traditional solvers currently in scipy. It is worth noting that the current solvers, except for lsmr(), date from as far back as the 1960’s. Maybe it is time to let in some newer technology.
COMMENT
I have thoroughly tested these routines. But I certainly hope this group will try interesting problems with it and push its limits. Please forward me any problems with unexpected results.
I look forward to your feedback.
Ron Jones
rejones7 at msn.com<mailto:Rejones7 at msn.com>
http://www.rejones7.net/


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20201012/1375a0af/attachment-0001.html>