New keyword name for linalg.lu
Dear all, In [1], I have proposed a new keyword for scipy.linalg.lu function, with prospective name "p_indices" as a bool switch. This will have the effect of having the permutation matrix P in A = P L U decomposition returned as a vector of indices such as [3, 1, 2, 0] meaning first and fourth row is swapped. With this you can still compute A through L[P, :] @ U without allocating a very costly 2D matrix of 1s and 0s. In tall arrays, say 400x10, the resulting 2D array P is 400x400 and in the vector form this reduces to 400x1. Hence, the switch to control this behavior is done by a new keyword and I don't know if "p_indices" is a nice one (doesn't feel like it). I considered, return_indices, perm_as_1d, p_as_1d, p_as_vec and bunch more but either too many underscores or too verbose. Thus, I'd like to ask for either validation or new ideas. We should maybe ping the NumPy audience for this too as this has no counterpart in NumPy. Thanks, ilhan [1] : https://github.com/scipy/scipy/pull/18358
FWIW, I'd suggest to reconsider adding the keyword at all. Maybe return a permutation vector (ideally, a LAPACK convention) and add a helper function to convert it to a permutation matrix? On Mon, Apr 24, 2023 at 12:41 PM Ilhan Polat <ilhanpolat@gmail.com> wrote:
Dear all,
In [1], I have proposed a new keyword for scipy.linalg.lu function, with prospective name "p_indices" as a bool switch.
This will have the effect of having the permutation matrix P in A = P L U decomposition returned as a vector of indices such as [3, 1, 2, 0] meaning first and fourth row is swapped. With this you can still compute A through L[P, :] @ U without allocating a very costly 2D matrix of 1s and 0s. In tall arrays, say 400x10, the resulting 2D array P is 400x400 and in the vector form this reduces to 400x1.
Hence, the switch to control this behavior is done by a new keyword and I don't know if "p_indices" is a nice one (doesn't feel like it). I considered, return_indices, perm_as_1d, p_as_1d, p_as_vec and bunch more but either too many underscores or too verbose.
Thus, I'd like to ask for either validation or new ideas. We should maybe ping the NumPy audience for this too as this has no counterpart in NumPy.
Thanks, ilhan
[1] : https://github.com/scipy/scipy/pull/18358 _______________________________________________ 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
That would break a lot of code, I'm afraid. We can deprecate it possibly, but I'd say that would take a lot of time. The helper function is actually -> np.eye(n)[P, :] that's probably why mathematica doesn't even bother with the matrix. On Mon, Apr 24, 2023 at 3:51 PM Evgeni Burovski <evgeny.burovskiy@gmail.com> wrote:
FWIW, I'd suggest to reconsider adding the keyword at all. Maybe return a permutation vector (ideally, a LAPACK convention) and add a helper function to convert it to a permutation matrix?
On Mon, Apr 24, 2023 at 12:41 PM Ilhan Polat <ilhanpolat@gmail.com> wrote:
Dear all,
In [1], I have proposed a new keyword for scipy.linalg.lu function,
with prospective name "p_indices" as a bool switch.
This will have the effect of having the permutation matrix P in A = P L
U decomposition returned as a vector of indices such as [3, 1, 2, 0] meaning first and fourth row is swapped. With this you can still compute A through L[P, :] @ U without allocating a very costly 2D matrix of 1s and 0s. In tall arrays, say 400x10, the resulting 2D array P is 400x400 and in the vector form this reduces to 400x1.
Hence, the switch to control this behavior is done by a new keyword and
I don't know if "p_indices" is a nice one (doesn't feel like it). I considered, return_indices, perm_as_1d, p_as_1d, p_as_vec and bunch more but either too many underscores or too verbose.
Thus, I'd like to ask for either validation or new ideas. We should
maybe ping the NumPy audience for this too as this has no counterpart in NumPy.
Thanks, ilhan
[1] : https://github.com/scipy/scipy/pull/18358 _______________________________________________ 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
_______________________________________________ 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: ilhanpolat@gmail.com
On Mon, Apr 24, 2023 at 9:52 AM Evgeni Burovski <evgeny.burovskiy@gmail.com> wrote:
FWIW, I'd suggest to reconsider adding the keyword at all. Maybe return a permutation vector (ideally, a LAPACK convention) and add a helper function to convert it to a permutation matrix?
To be clear, `lu()` currently returns the full permutation matrix. The new proposed functionality is to control changing the format of this returned object, not to control whether or not something gets returned. -- Robert Kern
Potentially stupid question but with the new keyword how does this differ from lu_factor?
On Mon, Apr 24, 2023 at 3:23 PM Jake Bowhay <jb9.bowhay@gmail.com> wrote:
Potentially stupid question but with the new keyword how does this differ from lu_factor?
The primary purpose of `lu_factor` is to provide a packed representation of the decomposition that can be passed essentially opaquely to `lu_solve` so it can do the solve efficiently. While one can disentangle them into their P, L, and U with other functions, `lu_factor` keeps things in exactly the format that LAPACK needs to do `lu_solve` efficiently. Getting the P, L, and U separately is the main purpose of `lu`. Sometimes, you want to do other things with those individually than just do efficient solves. We have the same pattern of division of labor with `cho_factor/cho_solve` and `cholesky`, for example. As for the pivot matrix format, it's usually at least as convenient to have the permutation indices as it is to have the full matrix, so that would still fit within `lu`'s purpose, IMO. -- Robert Kern
Another use case, typically, is to factor once and use it I on multiple right hand sides for linear equations. You might not yet have access to the right hand side but don't want to pay decomposition cost each time for Ax = b1 Ax = b2 And so on. Basically a quick hack to accommodate the use case with a low-level function instead of providing a cached operator, done in the very early times. In hindsight it's easy to criticize but it filled an important void for quite some time. On Tue, Apr 25, 2023, 02:08 Jake Bowhay <jb9.bowhay@gmail.com> wrote:
Thanks for the explanation, that makes sense to me _______________________________________________ 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: ilhanpolat@gmail.com
It would be nice to add a quick note to the docs explaining when/which you should use. Currently both state "Compute pivoted LU decomposition of a matrix." which while true isn't very helpful for a user trying to decide which function to pick!
Sorry for reviving this thread, but indeed "p_indices" sounding uglier the more I look at it. Just to summarize what it does to save you from the thread above; when P, L, U = scipy.linalg.lu(A) is run, currently, P is returning a full 2D array. If A is a tall array say, (25, 5) then P is necessarily (25, 25). And it is just a permutaiton matrix, a row shuffled np.eye(25). Instead, you can ask with this new keyword to return that shuffle pattern. as a 1D array and hence P becomes (25, ) array. And there are more than rare cases that only L and U is needed for other purposes. So we would like to make that matrix vs array switch with a keyword. I thought "p_indices = True/False" was a good idea but it became a bit disturbing. Other matlab following software offers a string selector with "outputForm='matrix'/'vector'" which forces users to remember the options hence I am voting for a TrueFalse flag. However I can't come up with the right keyword. Could you please offer some alternatives even just for inspiration? Best, ilhan On Tue, Apr 25, 2023 at 9:34 AM Jake Bowhay <jb9.bowhay@gmail.com> wrote:
It would be nice to add a quick note to the docs explaining when/which you should use. Currently both state "Compute pivoted LU decomposition of a matrix." which while true isn't very helpful for a user trying to decide which function to pick! _______________________________________________ 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: ilhanpolat@gmail.com
On Sun, May 21, 2023, at 17:08, Ilhan Polat wrote:
Sorry for reviving this thread, but indeed "p_indices" sounding uglier the more I look at it. Just to summarize what it does to save you from the thread above;
when
P, L, U = scipy.linalg.lu(A)
is run, currently, P is returning a full 2D array. If A is a tall array say, (25, 5) then P is necessarily (25, 25). And it is just a permutaiton matrix, a row shuffled np.eye(25). Instead, you can ask with this new keyword to return that shuffle pattern. as a 1D array and hence P becomes (25, ) array. And there are more than rare cases that only L and U is needed for other purposes.
So we would like to make that matrix vs array switch with a keyword. I thought "p_indices = True/False" was a good idea but it became a bit disturbing. Other matlab following software offers a string selector with "outputForm='matrix'/'vector'" which forces users to remember the options hence I am voting for a TrueFalse flag. However I can't come up with the right keyword.
Could you please offer some alternatives even just for inspiration?
How about reversing the meaning of the flag and naming it something like `p_full` or `full_p` with a `True` default? I also considered "squareform" instead of "full" as in `spatial.distance.squareform`, because there's a similar 1d <-> 2d square thing here, but the semantics are different enough that this name would probably be more confusing than helpful. András
Best, ilhan
On Tue, Apr 25, 2023 at 9:34 AM Jake Bowhay <jb9.bowhay@gmail.com> wrote:
It would be nice to add a quick note to the docs explaining when/which you should use. Currently both state "Compute pivoted LU decomposition of a matrix." which while true isn't very helpful for a user trying to decide which function to pick! _______________________________________________ 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: ilhanpolat@gmail.com
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: deak.andris@gmail.com
In the end I went with "matrix_p = True" to ring like "permute_l" which was I think what you were suggesting anyways. Hopefully some day in the future we can change this to False and be done with it. On Sun, May 21, 2023 at 11:31 PM Andras Deak <deak.andris@gmail.com> wrote:
On Sun, May 21, 2023, at 17:08, Ilhan Polat wrote:
Sorry for reviving this thread, but indeed "p_indices" sounding uglier the more I look at it. Just to summarize what it does to save you from the thread above;
when
P, L, U = scipy.linalg.lu(A)
is run, currently, P is returning a full 2D array. If A is a tall array say, (25, 5) then P is necessarily (25, 25). And it is just a permutaiton matrix, a row shuffled np.eye(25). Instead, you can ask with this new keyword to return that shuffle pattern. as a 1D array and hence P becomes (25, ) array. And there are more than rare cases that only L and U is needed for other purposes.
So we would like to make that matrix vs array switch with a keyword. I thought "p_indices = True/False" was a good idea but it became a bit disturbing. Other matlab following software offers a string selector with "outputForm='matrix'/'vector'" which forces users to remember the options hence I am voting for a TrueFalse flag. However I can't come up with the right keyword.
Could you please offer some alternatives even just for inspiration?
How about reversing the meaning of the flag and naming it something like `p_full` or `full_p` with a `True` default? I also considered "squareform" instead of "full" as in `spatial.distance.squareform`, because there's a similar 1d <-> 2d square thing here, but the semantics are different enough that this name would probably be more confusing than helpful.
András
Best, ilhan
On Tue, Apr 25, 2023 at 9:34 AM Jake Bowhay <jb9.bowhay@gmail.com>
wrote:
It would be nice to add a quick note to the docs explaining when/which you should use. Currently both state "Compute pivoted LU decomposition of a matrix." which while true isn't very helpful for a user trying to decide which function to pick! _______________________________________________ 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: ilhanpolat@gmail.com
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: deak.andris@gmail.com
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: ilhanpolat@gmail.com
Conceptually this seems related to sparse matrices so I wonder if there is some helpful terminology there. Your 1D option would return “column indices” for the non-zero entries in each row of the 2D option. Best, Todd
On 21 May 2023, at 16:08, Ilhan Polat <ilhanpolat@gmail.com> wrote:
[snip]
when
P, L, U = scipy.linalg.lu <http://scipy.linalg.lu/>(A)
is run, currently, P is returning a full 2D array. If A is a tall array say, (25, 5) then P is necessarily (25, 25). And it is just a permutaiton matrix, a row shuffled np.eye(25). Instead, you can ask with this new keyword to return that shuffle pattern. as a 1D array and hence P becomes (25, ) array. [snip] Could you please offer some alternatives even just for inspiration?
Best, ilhan
On Tue, Apr 25, 2023 at 9:34 AM Jake Bowhay <jb9.bowhay@gmail.com <mailto:jb9.bowhay@gmail.com>> wrote: It would be nice to add a quick note to the docs explaining when/which you should use. Currently both state "Compute pivoted LU decomposition of a matrix." which while true isn't very helpful for a user trying to decide which function to pick! _______________________________________________ SciPy-Dev mailing list -- scipy-dev@python.org <mailto:scipy-dev@python.org> To unsubscribe send an email to scipy-dev-leave@python.org <mailto:scipy-dev-leave@python.org> https://mail.python.org/mailman3/lists/scipy-dev.python.org/ <https://mail.python.org/mailman3/lists/scipy-dev.python.org/> Member address: ilhanpolat@gmail.com <mailto:ilhanpolat@gmail.com> _______________________________________________ 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: shoppert@baileywick.plus.com
Yes indeed that's something I checked but in sparse case there is no possibility of full array return due to the size constraints hence, say in SuperLU, they only have "perm_r" and "perm_c" integer 1D arrays anyways for row and column permutation indices. Also permutation index is a thing in combinatorics so I have p_as_vector = False/True return_indices = False/True so far that doesn't sound confusing and also not that ugly. On Mon, May 22, 2023 at 8:37 AM Todd Bailey <shoppert@baileywick.plus.com> wrote:
Conceptually this seems related to sparse matrices so I wonder if there is some helpful terminology there. Your 1D option would return “column indices” for the non-zero entries in each row of the 2D option.
Best, Todd
On 21 May 2023, at 16:08, Ilhan Polat <ilhanpolat@gmail.com> wrote:
[snip]
when
P, L, U = scipy.linalg.lu(A)
is run, currently, P is returning a full 2D array. If A is a tall array say, (25, 5) then P is necessarily (25, 25). And it is just a permutaiton matrix, a row shuffled np.eye(25). Instead, you can ask with this new keyword to return that shuffle pattern. as a 1D array and hence P becomes (25, ) array. [snip] Could you please offer some alternatives even just for inspiration?
Best, ilhan
On Tue, Apr 25, 2023 at 9:34 AM Jake Bowhay <jb9.bowhay@gmail.com> wrote:
It would be nice to add a quick note to the docs explaining when/which you should use. Currently both state "Compute pivoted LU decomposition of a matrix." which while true isn't very helpful for a user trying to decide which function to pick! _______________________________________________ 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: ilhanpolat@gmail.com
_______________________________________________ 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: shoppert@baileywick.plus.com
_______________________________________________ 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: ilhanpolat@gmail.com
That breaks practically all LU code out there so we have to keep the array return for quite a while before we attempt deprecating it. On Mon, May 22, 2023 at 3:33 PM Christian Lorentzen <lorentzen.ch@gmail.com> wrote:
Hi How about returning the most efficient one, probably a 1-dim ndarray and then providing good examples and/or utility functions for this returned object/ndarray? This way, the LU function does not get overloaded.
Best Christian
Am 22.05.2023 um 08:58 schrieb Ilhan Polat <ilhanpolat@gmail.com>:
Yes indeed that's something I checked but in sparse case there is no possibility of full array return due to the size constraints hence, say in SuperLU, they only have "perm_r" and "perm_c" integer 1D arrays anyways for row and column permutation indices. Also permutation index is a thing in combinatorics so I have
p_as_vector = False/True return_indices = False/True
so far that doesn't sound confusing and also not that ugly.
On Mon, May 22, 2023 at 8:37 AM Todd Bailey <shoppert@baileywick.plus.com> wrote:
Conceptually this seems related to sparse matrices so I wonder if there is some helpful terminology there. Your 1D option would return “column indices” for the non-zero entries in each row of the 2D option.
Best, Todd
On 21 May 2023, at 16:08, Ilhan Polat <ilhanpolat@gmail.com> wrote:
[snip]
when
P, L, U = scipy.linalg.lu(A)
is run, currently, P is returning a full 2D array. If A is a tall array say, (25, 5) then P is necessarily (25, 25). And it is just a permutaiton matrix, a row shuffled np.eye(25). Instead, you can ask with this new keyword to return that shuffle pattern. as a 1D array and hence P becomes (25, ) array. [snip] Could you please offer some alternatives even just for inspiration?
Best, ilhan
On Tue, Apr 25, 2023 at 9:34 AM Jake Bowhay <jb9.bowhay@gmail.com> wrote:
It would be nice to add a quick note to the docs explaining when/which you should use. Currently both state "Compute pivoted LU decomposition of a matrix." which while true isn't very helpful for a user trying to decide which function to pick! _______________________________________________ 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: ilhanpolat@gmail.com
_______________________________________________ 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: shoppert@baileywick.plus.com
_______________________________________________ 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: ilhanpolat@gmail.com
_______________________________________________ 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: lorentzen.ch@gmail.com
_______________________________________________ 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: ilhanpolat@gmail.com
participants (7)
-
Andras Deak -
Christian Lorentzen -
Evgeni Burovski -
Ilhan Polat -
Jake Bowhay -
Robert Kern -
Todd Bailey