scipy.linalg.lu considerations
Dear all, As I recently figured out (again) that SciPy has some custom Fortran code for generating LU decomposition and determinant [1], I started decommissioning those parts as we are slowly enabling nDarray support and some performance tweaks in the meantime [2]. I am about to post a PR for scipy.linalg.lu and it occured to me that we still have the legacy behavior that returns different number of outputs depending on the parameters (from the matlab inspired times). We painfully learned in the past not to do this anymore but we still have this behavior in parts of the API. I can definitely keep this behavior with no problem. Hence it is not a matter of possibility. But; As everybody knows, lu performs the decomposition A = P L U with P permutation matrix, lower and upper trapezoidal arrays. Now, forming a P as a complete and very sparsely populated array is a waste of memory if you are not interested and hence other software offers to return it as a 1D permutation array (SciPy doesn't but it is trivial to add this ala matlab [3] ) So proposition 1: Adding the output_form (or some other kwarg name) If the input is ndarray with shape (..., n, m) P becomes (..., min(m, n)) in the output_form='vector' case. Proposition 2: Going through the very annoying process of deprecating the mutating output arg number and return always p, l, u. Note that this is a central function in the number crunching community as there is no substitute in NumPy. I want to get some feedback from the folks about this to assess the weight of these propositions. Proposition 1 is basically just informing, I am +1 with the second one but definitely not looking forward to it :) Let me know, Best, ilhan [1]: https://github.com/scipy/scipy/blob/main/scipy/linalg/src/lu.f [2]: https://github.com/scipy/scipy/pull/18225 [3]: https://www.mathworks.com/help/matlab/ref/lu.html#d124e909668
Oh and just to add, if permute_l=True is given "P" would return empty but still 3 output arguments. On Thu, Apr 6, 2023 at 2:17 PM Ilhan Polat <ilhanpolat@gmail.com> wrote:
Dear all,
As I recently figured out (again) that SciPy has some custom Fortran code for generating LU decomposition and determinant [1], I started decommissioning those parts as we are slowly enabling nDarray support and some performance tweaks in the meantime [2].
I am about to post a PR for scipy.linalg.lu and it occured to me that we still have the legacy behavior that returns different number of outputs depending on the parameters (from the matlab inspired times). We painfully learned in the past not to do this anymore but we still have this behavior in parts of the API. I can definitely keep this behavior with no problem. Hence it is not a matter of possibility. But;
As everybody knows, lu performs the decomposition A = P L U with P permutation matrix, lower and upper trapezoidal arrays. Now, forming a P as a complete and very sparsely populated array is a waste of memory if you are not interested and hence other software offers to return it as a 1D permutation array (SciPy doesn't but it is trivial to add this ala matlab [3] )
So proposition 1: Adding the output_form (or some other kwarg name) If the input is ndarray with shape (..., n, m) P becomes (..., min(m, n)) in the output_form='vector' case.
Proposition 2: Going through the very annoying process of deprecating the mutating output arg number and return always p, l, u. Note that this is a central function in the number crunching community as there is no substitute in NumPy.
I want to get some feedback from the folks about this to assess the weight of these propositions. Proposition 1 is basically just informing, I am +1 with the second one but definitely not looking forward to it :)
Let me know,
Best, ilhan
[1]: https://github.com/scipy/scipy/blob/main/scipy/linalg/src/lu.f [2]: https://github.com/scipy/scipy/pull/18225 [3]: https://www.mathworks.com/help/matlab/ref/lu.html#d124e909668
Without looking into the details, how about using some form of what scipy stats does (or was doing?) with make_tuple_bunch --- basically a class which behaves like a tuple of a fixed number of arguments, and then P or a permutation matrix becomes an extra attribute? That said, why do we carry custom LU decomposition code at all, is this not a single LAPACK call (maybe plus a thin wrapper to preserve backcompat). On Thu, Apr 6, 2023 at 3:28 PM Ilhan Polat <ilhanpolat@gmail.com> wrote:
Oh and just to add, if permute_l=True is given "P" would return empty but still 3 output arguments.
On Thu, Apr 6, 2023 at 2:17 PM Ilhan Polat <ilhanpolat@gmail.com> wrote:
Dear all,
As I recently figured out (again) that SciPy has some custom Fortran code for generating LU decomposition and determinant [1], I started decommissioning those parts as we are slowly enabling nDarray support and some performance tweaks in the meantime [2].
I am about to post a PR for scipy.linalg.lu and it occured to me that we still have the legacy behavior that returns different number of outputs depending on the parameters (from the matlab inspired times). We painfully learned in the past not to do this anymore but we still have this behavior in parts of the API. I can definitely keep this behavior with no problem. Hence it is not a matter of possibility. But;
As everybody knows, lu performs the decomposition A = P L U with P permutation matrix, lower and upper trapezoidal arrays. Now, forming a P as a complete and very sparsely populated array is a waste of memory if you are not interested and hence other software offers to return it as a 1D permutation array (SciPy doesn't but it is trivial to add this ala matlab [3] )
So proposition 1: Adding the output_form (or some other kwarg name) If the input is ndarray with shape (..., n, m) P becomes (..., min(m, n)) in the output_form='vector' case.
Proposition 2: Going through the very annoying process of deprecating the mutating output arg number and return always p, l, u. Note that this is a central function in the number crunching community as there is no substitute in NumPy.
I want to get some feedback from the folks about this to assess the weight of these propositions. Proposition 1 is basically just informing, I am +1 with the second one but definitely not looking forward to it :)
Let me know,
Best, ilhan
[1]: https://github.com/scipy/scipy/blob/main/scipy/linalg/src/lu.f [2]: https://github.com/scipy/scipy/pull/18225 [3]: https://www.mathworks.com/help/matlab/ref/lu.html#d124e909668
_______________________________________________ SciPy-Dev mailing list -- scipy-dev@python.org To unsubscribe send an email to scipy-dev-leave@python.org https://mail.python.org/mailman3/lists/scipy-dev.python.org/ Member address: evgeny.burovskiy@gmail.com
Without looking into the details, how about using some form of what scipy stats does (or was doing?) with make_tuple_bunch --- basically a class which behaves like a tuple of a fixed number of arguments, and then P or a permutation matrix becomes an extra attribute?
Yes, that's probably how we are going to do it, if, we decide to do it.
That said, why do we carry custom LU decomposition code at all, is this not a single LAPACK call (maybe plus a thin wrapper to preserve backcompat).
That code is going away regardless. Det is done, with LU we can remove _flinalg.py stuff.
participants (2)
-
Evgeni Burovski -
Ilhan Polat