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