add more arguments to the function sfunction in linalg.ordqz

[freude]

In the current scipy version, one may pass a callable as the argument sort of the function linalg.ordqz. This callable is supposed to be a function taking two arguments - alpha and betha, which are thought to be enough to figure out the sorting procedure. I have found applications where sorting is performed is a way that it also requires information on the eigenvectors Q and Z. Therefore, I propose to spread the argument list of the callable making it to accept also Q and Z outputing from the function _qz - that function is normally called before sorting. If you find this feature useful - let me know and I will make a pull request. In the pull request I have changed the line 358 in the module _decomp_qz.py

from

select = sfunction(alpha, beta)

to

select = sfunction(alpha, beta, Q, Z)

[ilayn]

Can you give an example of such case? What kind of condition do you need such that you use QZ?

[freude]

Yes, here is an example of the problem. Each eigenvalue, E, can be computed as E=alpha/betha. In my case the eigenvalues are complex numbers and I need to sort them: if |E| < 1 - they go to the upper-left block, if |E| > 1 - lower-right block. So far, so good. The problem arises when |E| = 1 - classification of eigenvalues located on the unit circle in the complex plane. In this case, the decision can made using eigenvectors by computing so called group velocity of the wave packet in my example: Im[Z[:, j] * Z[:, j].T] > 0 - upper block, Im[Z[:, j] * Z[:, j].T] < 0 - lower block.

The modification I have mentioned can be made in a way that it does not affect the code requiring only alpha and betha by providing default values for Q and Z.

[ilayn]

This still doesn't sort unit sized eigs. Comparison would still be the same for all |E| =1 and to in my opinion I think this can be done much quicker after obtaining the QZ decomposition and working with the subset of blocks.

[freude]

Thank you for the response.

This still doesn't sort unit sized eigs.

How does this statement come? In my code it works perfectly. It sorts those eigenvalues in the way, that some of them with |E|=1 appear in the upper block, others - in the lower, and it is matters which one exactly goes where. It simply means that if we can not sort them by absolute value we can still sort them by angle, eigenvector properties or whatever.

Why do you think that only eigenvalues but not eigenvectors can be used to rearrange Shur forms - isn't it a lack of generality? I just don't want to patch scipy and add it into my package every time I need to distribute my software.

[freude]

In LAPACK sorting is done by the function ?TGSEN. It accepts a logical argument select which is fully user defined and library itself doesn't impose restrictions that it can be computed only on the basis of eigenvalues. However the complete information related to the problem is a set of eigenvalues and eigenvectors.

[ilayn]

So if you have two blocks both satisfying Im[Z[:, j] * Z[:, j].T] > 0 do you check their positiveness magnitude to sort them further? Or you just want them to appear at the top without any subsorting among |E|=1. That was my previous remark maybe worded a bit cryptic.

[freude]

First I check the magnitude of the eigenvalues, and only if the magnitude equals one I need further information from eigenvectors to define which block do they belong to. Let say I have 4 eigenvalues with |E|=1 - they always come in pairs - so I need to define 2 of them which go to the upper block and other 2 which go to the lower block.

Here is example of spectrum of absolute values: |E| = [ 0.1 0.5 1 1 1 1 2 10] - they are already sorted relative to the absolute value, but the order of those having |E|=1 may be wrong.

[freude]

So what is the final decision on the issue?

[rgommers]

@ilayn any more thoughts?

@freude maybe share your code (as a link, gist or open a PR) so it's easier to discuss about?