Skip to content

Comments

Replace Brandts sorting with LAPACK DTRSEN, stabilize Gram-Schmidt#64

Open
Marius1311 wants to merge 1 commit intomainfrom
feature/scipy-sorted-schur
Open

Replace Brandts sorting with LAPACK DTRSEN, stabilize Gram-Schmidt#64
Marius1311 wants to merge 1 commit intomainfrom
feature/scipy-sorted-schur

Conversation

@Marius1311
Copy link
Collaborator

  • Replace 624-line Python Brandts sorting (_sort_real_schur.py) with a single LAPACK DTRSEN call in sorted_brandts_schur(). Same Bai-Demmel algorithm, compiled Fortran instead of Python loops. Fixes subspace angle failures for n >= 100.

  • Replace single-pass Modified Gram-Schmidt in _gram_schmidt_mod() with Householder QR (np.linalg.qr, LAPACK DGEQRF). Single-pass MGS lost orthogonality for ill-conditioned inputs (kappa up to 1e16). QR is backward stable regardless of conditioning. Same O(nm^2) complexity, faster in practice (compiled LAPACK vs Python loops).

  • Replace fragile constant-vector detection (element-wise constancy check with tight tolerance) with cosine similarity against sqrt(eta). The old check failed when Schur vectors were only approximately constant due to numerical noise from different eigenvalue orderings.

  • Recompute R in _do_schur() after re-orthogonalization to maintain the Schur relation in the new basis.

  • Update tests accordingly.

- Replace 624-line Python Brandts sorting (_sort_real_schur.py) with a
  single LAPACK DTRSEN call in sorted_brandts_schur(). Same Bai-Demmel
  algorithm, compiled Fortran instead of Python loops. Fixes subspace
  angle failures for n >= 100.

- Replace single-pass Modified Gram-Schmidt in _gram_schmidt_mod() with
  Householder QR (np.linalg.qr, LAPACK DGEQRF). Single-pass MGS lost
  orthogonality for ill-conditioned inputs (kappa up to 1e16). QR is
  backward stable regardless of conditioning. Same O(nm^2) complexity,
  faster in practice (compiled LAPACK vs Python loops).

- Replace fragile constant-vector detection (element-wise constancy
  check with tight tolerance) with cosine similarity against sqrt(eta).
  The old check failed when Schur vectors were only approximately
  constant due to numerical noise from different eigenvalue orderings.

- Recompute R in _do_schur() after re-orthogonalization to maintain the
  Schur relation in the new basis.

- Update tests accordingly.
@Marius1311 Marius1311 force-pushed the feature/scipy-sorted-schur branch from 21ace14 to bd49915 Compare February 13, 2026 19:34
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant