BUG: sorting by imaginary part does not work for real Schur form in scipy.linalg.schur

[kaathaarinaa]

Describe your issue.

For a 3x3 matrix A I want to find its sorted real Schur form T ((quasi-)upper-triangular matrix) and the corresponding unitary Schur transformation matrix Z. I want the sorting to be done according to a rule defined in the callable kwarg sort such that the complex conjugate eigenvalue pair (if existent) will be placed at the top left of the Schur form. Therefore I use:

sorting_criterion = lambda x: np.invert(np.isclose(x.imag, 0, rtol=1e-5))

This returns True if the eigenvalue has (a reasonable) complex component.

Schur decomposition can be constructed either in real or in complex form. However, if constructed using scipy.linalg.schur in real form, sorting by whether or not the eigenvalue was complex does not work. It will return sdim = 0 which means the sorting criterion did not apply to any of the eigenvalues, even though they are initially complex.
LAPACK gees() (which scipy uses (?)) suggests though, that it can order the eigenvalues on the diagonal of the real Schur form so that selected eigenvalues are at the top left. Liu et al. 2018 use the LAPACK dgees() to get the real Schur form sorted in a way so that the real part of the complex conjugate eigenvalue pair appears on the top left of the Schur form.

Is the fact, that in scipy.linalg.schur the real Schur form cannot be sorted by the sorting_criterion by design or by accident?
Or am I missing something and this can already be done in scipy.linalg.schur too?

Example code with my sorting criterion below.

Reproducing Code Example

import numpy as np
from scipy import linalg as la

# A matrix to decompose
A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]])
# Decompose
sorting_criterion = lambda x: np.invert(np.isclose(x.imag, 0, rtol=1e-5))   # True if eigenvalue is imaginary
S, Q_star, sdim = la.schur(A, output='complex', sort=sorting_criterion)     # Complex decomposition
print(f'This is the sorting of the diagonal I want: {S.diagonal()}')
S_r, Q_star_r, sdim_r = la.schur(A, output='real', sort=sorting_criterion)  # Real decomposition
print(f'When using same sorting criterion in real Schur decomposition it yields: {S_r.diagonal()}')

Error message

no error message

SciPy/NumPy/Python version information

1.7.1 1.20.3 sys.version_info(major=3, minor=8, micro=2, releaselevel='final', serial=0)

[hlgsen]

I had a similar problem, but when sorting for absolute values of complex eigenvalues using sort='iuc'. This code produced my problem:

import numpy as np
from scipy.linalg import schur
A = np.array([[1, -0.5, 0], [2, 1, 0], [0, 0, 0.5]])
T, Z, sdim = schur(A, output='real', sort='iuc')
print(f"eigvals: {np.linalg.eigvals(A)}")
print(f"sdim suggests {sdim} eigenvalues with |lambda| <= 1, correct sdim would be 1")
print(f"T = {T}, notice lambda=0.5 is not on top left")

I found this related issue. Passing a function with two arguments

def cond(x_re, x_im):
    return x_re ** 2 + x_im ** 2 < 1

using sort=cond was a workaround to my problem. In your case, maybe something like

def sorting_criterion(x_re, x_im):
    return  np.invert(np.isclose(x_im, 0, rtol=1e-5))

would work?