[Numpy-discussion] Fwd: Backslash operator A\b and np/sp.linalg.solve

Mon Jan 9 05:33:56 EST 2017

On Mon, Jan 9, 2017 at 4:17 AM, Ilhan Polat <ilhanpolat at gmail.com> wrote:

>
> Hi everyone,
>
> I was stalking the deprecating the numpy.matrix discussion on the other
> thread and I wondered maybe the mailing list is a better place for the
> discussion about something I've been meaning to ask the dev members. I
> thought mailing lists are something we dumped using together with ICQ and
> geocities stuff but apparently not :-)
>
> Anyways, first thing is first: I have been in need of the ill-conditioned
> warning behavior of matlab (and possibly other software suites) for my own
> work. So I looked around in the numpy issues and found
> https://github.com/numpy/numpy/issues/3755 some time ago. Then I've
> learned from @rkern that there were C translations involved in the numpy
> source and frankly I couldn't even find the entry point of how the project
> is structured so I've switched to SciPy side where things are a bit more
> convenient. Next to teaching me more about f2py machinery, I have noticed
> that the linear algebra module is a bit less competitive than the usual
> level of scipy though it is definitely a personal opinion.
>
> So in order to get the ill-conditioning (or at least the condition number)
> I've wrapped up a PR using the expert routines of LAPACK (which is I think
> ready to merge) but still it is far from the contemporary software
> convenience that you generally get.
>
> https://github.com/scipy/scipy/pull/6775
>
> The "assume_a" keyword introduced here is hopefully modular enough that
> should there be any need for more structures we can simply keep adding to
> the list without any backwards compatibility. It will be at least offering
> more options than what we have currently. The part that I would like to
> discuss requires a bit of intro so please bear with me. Let me copy/paste
> the part from the old PR:
>
>   Around many places online, we can witness the rant about numpy/scipy not
> letting the users know about the conditioning for example Mike Croucher's
> blog <http://www.walkingrandomly.com/?p=5092> and numpy/numpy#3755
> <https://github.com/numpy/numpy/issues/3755>
>
>   Since we don't have any central backslash function that optimizes
> depending on the input data, should we create a function, let's name it
> with the matlab equivalent for the time being linsolve such that it
> automatically calls for the right solver? This way, we can avoid writing
> new functions for each LAPACK driver . As a reference here is a SO thread
> <http://stackoverflow.com/questions/18553210/how-to-implement-matlabs-mldivide-a-k-a-the-backslash-operator>
> that summarizes the linsolve
> <http://nl.mathworks.com/help/matlab/ref/linsolve.html> functionality.
>

Note that you're proposing a new scipy feature (right?) on the numpy
list....

This sounds like a good idea to me. As a former heavy Matlab user I
remember a lot of things to dislike, but "\" behavior was quite nice.

> I'm sure you are aware, but just for completeness, the linear equation
> solvers are often built around the concept of polyalgorithm which is a
> fancy way of saying that the array is tested consecutively for certain
> structures and the checks are ordered in such a way that the simpler
> structure is tested the sooner. E.g. first check for diagonal matrix, then
> for upper/lower triangular then permuted triangular then symmetrical and so
> on. Here is also another example from AdvanPix
> http://www.advanpix.com/2016/10/07/architecture-of-linear-systems-solver/
>
> Now, according to what I have coded and optimized as much as I can, a pure
> Python is not acceptable as an overhead during these checks. It would
> definitely be a noticeable slowdown if this was in place in the existing
> linalg.solve however I think this is certainly doable in the low-level
> C/FORTRAN level.
>

How much is a noticeable slowdown? Note that we still have the current
interfaces available for users that know what they need, so a nice
convenience function that is say 5-10% slower would not be the end of the
world.

Ralf

> CPython is certainly faster but I think only a straight C/FORTRAN
> implementation would cut it. Note that we only need the discovery of the
> structure then we can pass to the dedicated solver as is. Hence I'm not
> saying that we should recode the existing solve functionality. We already
> have the branching in place to ?GE/SY/HE/POSVX routines.
>
> -------
>
> The second issue about the current linalg.solve function is when trying to
> solve for right inverse e.g. xA = b. Again with some copy/paste: The right
> inversion is currently a bit annoying, that is to say if we would like to
> compute, say, BA^{-1}, then the user has to explicitly transpose the
> explicitly transposed equation to avoid using an explicit inv(whose use
> should be discouraged anyways)
> x = scipy.linalg.solve(A.T, B.T).T.
>
> Since expert drivers come with a trans switch that can internally handle
> whether to solve the transposed or the regular equation, these routines
> avoid the A.T off-the-shelf. I am wondering what might be the best way to
> add a "r_inv" keyword such that the B.T is also handled at the FORTRAN
> level instead such that the user can simply write "solve(A,B, r_inv=True)".
> Because we don't have a backslash operation we could at least provide this
> much as convenience I guess.
>
> I would love to have go at it but I'm definitely not competent enough in
> C/FORTRAN at the production level so I was wondering whether I could get
> some help about this. Anyways, I hope I could make my point with a rather
> lengthy post. Please let me know if this is a plausible feature
>
> ilhan
>
> PS: In case gmail links won't be parsed, here are the inline links
>
> MC blog: http://www.walkingrandomly.com/?p=5092
> SO thread : http://stackoverflow.com/questions/18553210/how-to-implement
> -matlabs-mldivide-a-k-a-the-backslash-operator
> linsolve/mldivide page : http://nl.mathworks.com/help/m
> atlab/ref/mldivide.html
>
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20170109/84c2fd2a/attachment.html>