Mailman 3 October 2020 - SciPy-Dev

Resubmission of autosol.py
by rondall jones Oct. 13, 2020

Oct. 13, 2020

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(a)msn.com<mailto:Rejones7@msn.com> http://www.rejones7.net/

3 7

PR Idea: Allowing multiple axis arguments (axis = tuple of ints) in the stats package?
by Thomas Hodson Oct. 6, 2020

Oct. 6, 2020

Hello! I use scipy.stats.sem a lot and I would love if it was able to take multiple axis arguments as many numpy functions can. Looking at the source code sem and other functions in scipy.stats are already implemented mostly in terms of numpy functions so it seems like it would only require changing some parts of the logic. Taking stats.sem as an example I think it would require changing n = a.shape[axis] to something like n = product(a.shape[axis]) For masked arrays, a.count(axis) is used which already works with multiple axes. Before I start work on a PR I wanted to ask if there is some reason that this change would be considered a bad idea. Otherwise, if I write a decent PR with tests, benchmarks and documentation that updates all the functions in scipy.stats that might reasonably take multiple axis arguments, is it likely to be accepted? Thanks, Tom

2 1

Re: [SciPy-Dev] Proposal for a new ODE Solver
by Laurin Steidle Oct. 5, 2020

Oct. 5, 2020

Hi all, I've proposed the implementation of the modified-Patanka-Runge-Kutta (MPRK) ODE solver. (original message sent at Thu, 27 Aug 2020) I've continued with the development but now I am stuck at a point where decision needs to be made. I've presented more detailed in the github issue <https://github.com/scipy/scipy/issues/12766> but to summarize it quickly: MPRK requires the form of the ODE to be |dy / dt = p(t, y) + d(t, y) where (p/d) is (positive/negative) | instead of the scipy default: || |dy / dt = f(t, y) ||| This opens up three options on how to proceed from here (that I can think of): 1. Allow the new format to be used in solve_ivp. (creates two disjoint ODE solver sets in solve_ivp) 2. Create a new "solve_ivp"-like function. (creates a naming conflict as both methods solve ivp's) 3. Abandon the new solver I personally don't like either of these options very much but I am inclined to the first. Reasons for that are outlined in the github issue. Any thoughts on this? Best, Laurin

1 0

Re: [SciPy-Dev] Discontinued Rackspace Open Compute Discount
by Robert Lucente Oct. 5, 2020

Oct. 5, 2020

Several CI services are free for open source projects although there are throughput limits https://about.gitlab.com https://github.com/features/actions https://travis-ci.org https://www.appveyor.com -----Original Message----- >From: Olivier Grisel <olivier.grisel(a)ensta.org> >Sent: Feb 8, 2020 11:30 AM >To: SciPy Developers List <scipy-dev(a)python.org> >Subject: [SciPy-Dev] Discontinued Rackspace Open Compute Discount > >Hi all, > >I received an email end of January that I just discovered today >telling me that the Rackspace Open Compute Discount grant we use to >host http://wheels.scipy.org for nightly builds and release management >of our community has been discontinued on December 31st 2019. > >I would like to take down the Rackspace hosted cloudfiles as soon as >possible but we need to find a replacement as this will break all the >CIs that rely on nightly wheels and also all the wheel release >infrastructure for several project in the ecosystem. > >Apparently the fees for January and the beginning of February already >sum to $936.51. This is a large cost because there were a bunch of >unused left-over block volumes that I already deleted. > >To host wheels for nightly builds and release automation, a possible >alternatives include using the github release pages but they are not >really designed for this kind of workflow. > >Alternatively there is a PyPI repo hosting feature in Azure Pipelines >in public preview but I am not sure how it works, if it's free for >open source projects and if we could use it as a replacement for the >Rackspace cloud storage we currently use to host >http://wheels.scipy.org . > >https://docs.microsoft.com/en-us/azure/devops/pipelines/artifacts/pypi?view… > >Important: if you have uploaded any important files to one of the >cloud file storage / blob containers on Rackspace, please download a >copy as soon as possible. > >But please do not download things you do not really want to keep to >avoid adding any unnecessary bandwidth costs. If you have file write >permission on one of those containers, please delete the files you own >and do not need anymore. > >Feel free to join this room to discuss the matter interactively: >https://gitter.im/scikit-learn/nightly-build-hosting > >Feel free to reply to this email on the scipy-dev mailing list if you >have suggestions for setting up a free replacement. > >-- >Olivier >http://twitter.com/ogrisel - http://github.com/ogrisel >_______________________________________________ >SciPy-Dev mailing list >SciPy-Dev(a)python.org >https://mail.python.org/mailman/listinfo/scipy-dev

12 25

Enhance stats.shapiro() to take n-dimension input and nan_policy
by Po-Wen Kao Oct. 4, 2020

Oct. 4, 2020

Hi scipy-dev, Recently I am working on data analysis combining Pandas and scipy. I realize some API in stats module takes only 1d array and can't omit np.nan value. So I create a PR (https://github.com/scipy/scipy/pull/12916) on the function that I use the most. If the community finds this feasible, I would be happy to help also upgrade other APIs in the future. :) In case you find this PR interesting, feel free to be the reviewer of the change. Thanks Best Regards, Po-Wen Kao

1 0

Re-running Kmeans on a dataset in Python
by Stephen Malcolm Oct. 4, 2020

Oct. 4, 2020

Hello, I've written some code to run Kmeans on a data set (please see below). And I've plotted the results, with my two clusters/ centroids. However, I've to re-run Kmeans several times and pull up different plots (showing the different centroid positions). Can someone point me in the right direction how to write this extra code to perform this task? Then I've to conclude if Kmeans is stable. I believe this is the lowest sum of squared errors? Thanking you in advance. pandas used to read dataset and return the data #numpy and matplotlib to represent and visualize the data #sklearn to implement kmeans algorithm import pandas as pd import numpy as np import matplotlib.pyplot as plt from sklearn.cluster import KMeans from sklearn.preprocessing import StandardScaler #import the data data = pd.read_csv('file.csv') #extract values x=data['V1'] y=data['V2'] V1_V2 = np.column_stack ((V1, V2)) km_res = KMeans (n_clusters= 2).fit(V1_V2) y_kmeans = km_res.predict(V1_V2) plt.scatter(V1, V2, c=y_kmeans, cmap='viridis', s = 50, alpha = 0.5) plt.xlabel('V1') plt.ylabel('V2') plt.title('Visualization of raw data'); clusters = km_res.cluster_centers_ plt.scatter(clusters[:,0], clusters[:,1], c='blue', s=150)

1 0

Issues Inserting Graphical Overlay Using Matplotlib Patches
by Stephen Malcolm Oct. 1, 2020

Oct. 1, 2020

Hello All, I'm having some trouble adding a graphical overlay i.e. an ellipse onto my plot. I wish to do this, as I need to explain/ portray the mean, standard deviation and outliers. And hence evaluate the suitability of the dataset. Could you please let me know what code I'm missing/ or need to add, in order to insert this ellipse? I have no trouble plotting the data points and the mean using this code, however, the ellipse (width and height/ standard deviation) doesn't appear. I have no errors, instead, I'm getting a separate graph (without data points or ellipse) below the plotted one. Please find my code below: #pandas used to read dataset and return the data #numpy and matplotlib to represent and visualize the data #sklearn to implement kmeans algorithm import pandas as pd import numpy as np import matplotlib.pyplot as plt from sklearn.cluster import KMeans #import the data data = pd.read_csv('banknotes.csv') #extract values x=data['V1'] y=data['V2'] #print range to determine normalization print ("X_max : ",x.max()) print ("X_min : ",x.min()) print ("Y_max : ",y.max()) print ("Y_min : ",y.min()) #normalize values mean_x=x.mean() mean_y=y.mean() max_x=x.max() max_y=y.max() min_x=x.min() min_y=y.min() for i in range(0,x.size): x[i] = (x[i] - mean_x) / (max_x - min_x) for i in range(0,y.size): y[i] = (y[i] - mean_y) / (max_y - min_y) #statistical analyis using mean and standard deviation import matplotlib.patches as patches mean = np.mean(data, 0) std_dev = np.std(data, 0) ellipse = patches.Ellipse([mean[0], mean [1]], std_dev[0]*2, std_dev[1]*2, alpha=0.25) plt.xlabel('V1') plt.ylabel('V2') plt.title('Visualization of raw data'); plt.scatter(data.iloc[:, 0], data.iloc[:, 1]) plt.scatter(mean[0],mean[1]) plt.figure(figsize=(6, 6)) fig,graph = plt.subplots() graph.add_patch(ellipse) Kind Regards, Stephen

3 5