From 30b9f56a07b94ee8c1ddcef9d195033f6528acc5 Mon Sep 17 00:00:00 2001 From: Oleksandr Date: Tue, 30 Jun 2026 06:36:34 +0200 Subject: [PATCH] Fix infinite loop in complex non-symmetric EVD on defective matrices The non-symmetric Hessenberg-QR eigenvalue loops for Complex and Complex32 had no iteration cap, so they could loop forever on defective matrices with repeated eigenvalues (e.g. the 3x3 from #11). The real-valued Double/Single versions already guard this with 'iter >= 30*order -> NonConvergenceException'; this ports the same guard into the four complex-valued spots (UserEvd and the managed provider, for both Complex and Complex32), and adds regression tests covering each path. Fixes #11. --- .../Complex/Factorization/EvdTests.cs | 21 +++++++++++++++++++ .../Complex/Factorization/UserEvdTests.cs | 20 ++++++++++++++++++ .../Complex32/Factorization/EvdTests.cs | 21 +++++++++++++++++++ .../Complex32/Factorization/UserEvdTests.cs | 20 ++++++++++++++++++ .../Complex/Factorization/UserEvd.cs | 4 ++++ .../Complex32/Factorization/UserEvd.cs | 4 ++++ .../ManagedLinearAlgebraProvider.Complex.cs | 4 ++++ .../ManagedLinearAlgebraProvider.Complex32.cs | 4 ++++ 8 files changed, 98 insertions(+) diff --git a/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs b/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs index 59a6a70c0..48f4bd335 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/EvdTests.cs @@ -259,5 +259,26 @@ public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven([Valu AssertHelpers.AlmostEqual(ACopy, A, 14); AssertHelpers.AlmostEqual(BCopy, B, 14); } + + /// + /// Factorizing a defective matrix with a repeated eigenvalue throws + /// NonConvergenceException instead of looping forever (see issue #11). + /// This exercises the managed provider path used by DenseMatrix.Evd(). + /// + [Test] + public void DefectiveMatrixThrowsNonConvergenceException() + { + // Lower-triangular with a single eigenvalue (0.9275488) of algebraic + // multiplicity 3 - a defective (non-diagonalizable) matrix that the + // non-symmetric complex QR iteration cannot converge on. + var matrix = Matrix.Build.DenseOfArray(new[,] + { + { new Complex(0.9275488, 0), Complex.Zero, Complex.Zero }, + { new Complex(-0.9275488, 0), new Complex(0.9275488, 0), Complex.Zero }, + { Complex.Zero, new Complex(-0.9275488, 0), new Complex(0.9275488, 0) } + }); + + Assert.Throws(() => matrix.Evd()); + } } } diff --git a/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs b/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs index 5a8dadc05..41f48a862 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Complex/Factorization/UserEvdTests.cs @@ -314,5 +314,25 @@ public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven([Valu AssertHelpers.AlmostEqual(ACopy, A, 14); AssertHelpers.AlmostEqual(BCopy, B, 14); } + + /// + /// Factorizing a defective matrix with a repeated eigenvalue throws + /// NonConvergenceException instead of looping forever (see issue #11). + /// + [Test] + public void DefectiveMatrixThrowsNonConvergenceException() + { + // Lower-triangular with a single eigenvalue (0.9275488) of algebraic + // multiplicity 3 - a defective (non-diagonalizable) matrix that the + // non-symmetric complex QR iteration cannot converge on. + var matrixA = new UserDefinedMatrix(new[,] + { + { new Complex(0.9275488, 0), Complex.Zero, Complex.Zero }, + { new Complex(-0.9275488, 0), new Complex(0.9275488, 0), Complex.Zero }, + { Complex.Zero, new Complex(-0.9275488, 0), new Complex(0.9275488, 0) } + }); + + Assert.Throws(() => matrixA.Evd()); + } } } diff --git a/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs b/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs index ff1cf8efd..be4747dbb 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/EvdTests.cs @@ -257,5 +257,26 @@ public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven([Valu AssertHelpers.AlmostEqual(ACopy, A, 14); AssertHelpers.AlmostEqual(BCopy, B, 14); } + + /// + /// Factorizing a defective matrix with a repeated eigenvalue throws + /// NonConvergenceException instead of looping forever (see issue #11). + /// This exercises the managed provider path used by DenseMatrix.Evd(). + /// + [Test] + public void DefectiveMatrixThrowsNonConvergenceException() + { + // Lower-triangular with a single eigenvalue (0.9275488) of algebraic + // multiplicity 3 - a defective (non-diagonalizable) matrix that the + // non-symmetric complex QR iteration cannot converge on. + var matrix = Matrix.Build.DenseOfArray(new[,] + { + { new Complex32(0.9275488f, 0f), Complex32.Zero, Complex32.Zero }, + { new Complex32(-0.9275488f, 0f), new Complex32(0.9275488f, 0f), Complex32.Zero }, + { Complex32.Zero, new Complex32(-0.9275488f, 0f), new Complex32(0.9275488f, 0f) } + }); + + Assert.Throws(() => matrix.Evd()); + } } } diff --git a/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs b/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs index d66777007..dccac157e 100644 --- a/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs +++ b/src/Numerics.Tests/LinearAlgebraTests/Complex32/Factorization/UserEvdTests.cs @@ -312,5 +312,25 @@ public void CanSolveForRandomMatrixAndSymmetricMatrixWhenResultMatrixGiven([Valu AssertHelpers.AlmostEqual(ACopy, A, 14); AssertHelpers.AlmostEqual(BCopy, B, 14); } + + /// + /// Factorizing a defective matrix with a repeated eigenvalue throws + /// NonConvergenceException instead of looping forever (see issue #11). + /// + [Test] + public void DefectiveMatrixThrowsNonConvergenceException() + { + // Lower-triangular with a single eigenvalue (0.9275488) of algebraic + // multiplicity 3 - a defective (non-diagonalizable) matrix that the + // non-symmetric complex QR iteration cannot converge on. + var matrixA = new UserDefinedMatrix(new[,] + { + { new Complex32(0.9275488f, 0f), Complex32.Zero, Complex32.Zero }, + { new Complex32(-0.9275488f, 0f), new Complex32(0.9275488f, 0f), Complex32.Zero }, + { Complex32.Zero, new Complex32(-0.9275488f, 0f), new Complex32(0.9275488f, 0f) } + }); + + Assert.Throws(() => matrixA.Evd()); + } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs index b8409882b..f0c714387 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserEvd.cs @@ -662,6 +662,10 @@ static void NonsymmetricReduceHessenberToRealSchur(Matrix eigenVectors, exshift += s; iter++; + if (iter >= 30*order) + { + throw new NonConvergenceException(); + } // Reduce to triangle (rows) for (var i = l + 1; i <= n; i++) diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs index d5dbca343..6fbd1366d 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserEvd.cs @@ -666,6 +666,10 @@ static void NonsymmetricReduceHessenberToRealSchur(Matrix eigenVector exshift += s; iter++; + if (iter >= 30*order) + { + throw new NonConvergenceException(); + } // Reduce to triangle (rows) for (var i = l + 1; i <= n; i++) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index 97ca54148..db3bcb4ee 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -3055,6 +3055,10 @@ internal static void NonsymmetricReduceHessenberToRealSchur(Complex[] vectorV, C exshift += s; iter++; + if (iter >= 30*order) + { + throw new NonConvergenceException(); + } // Reduce to triangle (rows) for (var i = l + 1; i <= n; i++) diff --git a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index fbbc96f08..d5d0e3ddf 100644 --- a/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Providers/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -3057,6 +3057,10 @@ internal static void NonsymmetricReduceHessenberToRealSchur(Complex32[] vectorV, exshift += s; iter++; + if (iter >= 30*order) + { + throw new NonConvergenceException(); + } // Reduce to triangle (rows) for (var i = l + 1; i <= n; i++)