Mailman 3 August 2018 - SciPy-Dev

Review of PR 8293 - sampling of random variates
by Christoph Baumgarten Aug. 4, 2018

Aug. 4, 2018

PR 8293 introduces a method to sample random variates for more complex distributions to scipy.stats (such as hyperbolic distributions and the generalized inverse gaussian (see PR 8681)). I think the PR is in good shape, it would be great if it could be merged soon since i could then continue to work on adding new distributions. The PR has been open for 7 months now. If someone could continue the review, that would bei great. Thanks!

4 4

Backwards Compatibility for low level LAPACK routines
by Ilhan Polat Aug. 4, 2018

Aug. 4, 2018

Due their historical evolution, there are certain LAPACK wrappers that are not standardized. Some work with minimum lwork variables instead of their optimal values. Also these routines often return quite big arrays during the lwork queries, to demonstrate : import scipy.linalg as la la.lapack.dgeqrf(a=np.random.rand(400,400), lwork=-1) is a workspace size query (via lwork=-1). The current default size is "3*a.shape[0] + 1" hence 1201. However the optimal workspace size is 12800 on my machine. Therefore the mismatch is sometimes quite dramatic especially in some other routines. Notice also that to obtain this number the routine actually returns a 400-long array tau and requires the input matrix to be transferred back and forth. Moreover, they can't be handled via scipy.linalg.lapack._compute_lwork function. There are a few routines like this and I feel like they should be fixed and I'm willing to. However this means that their output signature is going to change which imply backwards compatibility breaks. I tried to see whether we could deprecate them with new wrappers with modified names, but to be honest, that would create too many duplicates. On the other hand I don't have a feeling of how much break this would mean out there in the wild. Is this break an acceptable one or not? (well, none is acceptable preferably, but in despair...) Any other alternatives, thoughts are most welcome. best, ilhan

4 7

Re: [SciPy-Dev] SciPy-Dev Digest, Vol 177, Issue 21
by Matt Newville July 31, 2018

July 31, 2018

Hi Andrew, Dickson, All, Sorry for barging in rather late to this discussion, and for responding to a digest, and if this appears to be too negative. > Date: Mon, 23 Jul 2018 23:36:37 -0400 > From: dickson41152 <dickson41152(a)protonmail.com> > To: Andrew Nelson <andyfaff(a)gmail.com> > Cc: scipy-dev <scipy-dev(a)python.org> > Subject: Re: [SciPy-Dev] Proposed Modifications to Nelder-Mead Simplex > Algorithm > Message-ID: > > <5xgSJzGhTqz3FgLpZlSQCqNBm-v5aytNC9kOBDasi_Mus8NNwoCKjPVEmrVqOpIn2B22oT1r0jOc8cNWJaBxRicbjjkAClPDRUwCRdPgAXs=@ > protonmail.com> > > Content-Type: text/plain; charset="utf-8" > > Thank you for your interest Andrew. > > > variable transformations to form unconstrained problems from bound > constrained problems > > I am considering writing a new transformation function to calculate the > equations on that web page that you linked to discussing variable > transformations. If it is as straight-forward as I think it will be to > write, then I may also do some comparisons between the "clipping" method > and the sin() / arcsin() / sqrt() method you highlighted. > > As an aside, I was told informally about the existence of that method you > referenced during a conversation I had last week, but I have been finding > it very difficult to find any written discussion (i.e. actual analytical > expressions) about that method via web searches before I made the > proposals. So thank you very much for linking that web page - it is very > useful information for my own personal work as well as the proposed Scipy > modifications! > > > user-specified objective function values for points in the > user-specified initial simplex > > I have indeed used a keyword in the manner you describe (and therefore > similar in style to the initial_simplex keyword) - it seemed to me to be a > "clean" way of implementing this modification and very straightforward to > implement since it is so similar to the initial_simplex handling. > > > evaluation budget count check > > I am thinking of writing checks which will effectively skip the execution > of the remainder of the solver algorithm if the budget has been met, before > finally sorting the current simplex of evaluations and returning the best. > My inclination is to avoid using StopIteration since it is an exception > i.e. it suggests some kind of problem has occurred, where as my proposed > idea is to handle evaluations up to the exact budget as a feature. > > Anyway I haven't yet made an attempt on this modification, so my thoughts > are tentative and may change - especially if it turns out that I am trying > to write spaghetti. > > dickson I think that adding bounds to Nelder-Mead Simplex is an interesting idea. Using the Minuit-style bounds - using smooth functions to transform infiniite "internal" values to finite "external" values - would definitely be an easy and stable way to do that. But, with all due respect, I also think it is the wrong approach ;). The main point of the lmfit package that Andrew linked to is that the fitting *algorithm* should not support bounds on variable values. I do understand that other scipy.optimize algorithms support bounds, and yes, I do think that is an unfortunate choice ;). Instead, the *parameters* should have bounds (and perhaps other properties), and the algorithm should work on these parameter objects. In this way one does not add bounds to Nelder-Mead, and then add bounds to Powell's method, and then add bounds to Levenberg-Marquardt, etc. Rather, one adds bounds to parameters, and has the methods use these parameters instead of scalar values, and thus one add bounds to *all* methods in a consistent and transferrable way. A few months ago Andrew suggested a class-based approach to implementing fitting algorithms. I think this idea has a lot of merit to it. I also believe that lmfit successfully demonstrates that having parameter objects as described briefly above does better encapsulate the desired properties (bounds, etc) *and* makes the code implementing the algorithms clearer and simpler. In fact, if I understand the original goals here to be: 1. Allowing bound constraints 2. Allowing the user to specify initial objective function values for any or all points in a user's initial simplex. 3. Checking the objective function evaluation count against the user-specified evaluation budget at each function evaluation . I would point out that lmfit has a Nelder-Mead solver that already does #1 and #3 and that (as with bounds), #3 may be best done not within the solver, but in outer code that calls the solver, so that this may be available in a consistent manner to multiple methods. If #2 means to memoize each step in the process to avoid repeated calculation, this does not currently exist within lmfit. Such a goal might really be best done in the solver method itself. It may be also be doable using a per-iteration callback (which lmfit does have) , though I don't know that anyone one has ever tried this. I apologize if this appears to be an outsider yelling "you're doing it all wrong". I rely heavily on scipy.optimize for my daily work, and would like to be able to make these tools better but I am over-committed and cannot contribute significantly to the scipy ecosystem myself other than try to support lmfit as best I can. Cheers, --Matt

2 1