Eigenvalues: Start ComplexQZ sweeps below two small subdiagonals libeigen/eigen!3430 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/Eigenvalues/ComplexQZ.h b/Eigen/src/Eigenvalues/ComplexQZ.h index f3a0ff0..d177617 100644 --- a/Eigen/src/Eigenvalues/ComplexQZ.h +++ b/Eigen/src/Eigenvalues/ComplexQZ.h
@@ -242,7 +242,9 @@ return numext::abs(x) <= tol; } - void do_QZ_step(Index p, Index q, unsigned int iter); + void do_QZ_step(Index f, Index l, unsigned int iter); + + Vec3 shiftColumn(Index p, Index l, unsigned int iter) const; JacobiRotation<Scalar> computeZk2(const Row2& b); @@ -455,7 +457,7 @@ if (z >= f) { push_down_zero_ST(z, l); } else { - do_QZ_step(f, m_n - l - 1, local_iter); + do_QZ_step(f, l, local_iter); local_iter++; m_global_iter++; } @@ -476,34 +478,57 @@ return J.transpose(); } +/** \internal First column of the shift polynomial of the block (S, T)(p:l, p:l), up to a nonzero factor. */ template <typename MatrixType_> -void ComplexQZ<MatrixType_>::do_QZ_step(Index p, Index q, unsigned int iter) { - // This is certainly not the most efficient way of doing this, - // but a readable one. +typename ComplexQZ<MatrixType_>::Vec3 ComplexQZ<MatrixType_>::shiftColumn(Index p, Index l, unsigned int iter) const { const auto a = [p, this](Index i, Index j) { return m_S(p + i - 1, p + j - 1); }; const auto b = [p, this](Index i, Index j) { return m_T(p + i - 1, p + j - 1); }; - const Index m = m_n - p - q; // Size of the inner block - Scalar x, y, z; + const Index m = l - p + 1; // Size of the inner block if (iter > 0 && iter % 10 == 0) { // Break stalled double-shift iterations with a single shift displaced by the trailing subdiagonals. const RealScalar displacement = numext::abs(a(m, m - 1) / b(m - 1, m - 1)) + numext::abs(a(m - 1, m - 2) / b(m - 2, m - 2)); const Scalar shift = a(m, m) / b(m, m) + displacement; // (S*T^-1 - shift*I)*e1 has only two nonzero entries. - x = a(1, 1) / b(1, 1) - shift; - y = a(2, 1) / b(1, 1); - z = Scalar(0); - } else { - Scalar W1 = a(m - 1, m - 1) / b(m - 1, m - 1) - a(1, 1) / b(1, 1), W2 = a(m, m) / b(m, m) - a(1, 1) / b(1, 1), - W3 = a(m, m - 1) / b(m - 1, m - 1); - - x = (W1 * W2 - a(m - 1, m) / b(m, m) * W3 + W3 * b(m - 1, m) / b(m, m) * a(1, 1) / b(1, 1)) * b(1, 1) / a(2, 1) + - a(1, 2) / b(2, 2) - a(1, 1) / b(1, 1) * b(1, 2) / b(2, 2); - y = (a(2, 2) / b(2, 2) - a(1, 1) / b(1, 1)) - a(2, 1) / b(1, 1) * b(1, 2) / b(2, 2) - W1 - W2 + - W3 * (b(m - 1, m) / b(m, m)); - z = a(3, 2) / b(2, 2); + return Vec3(a(1, 1) / b(1, 1) - shift, a(2, 1) / b(1, 1), Scalar(0)); } + const Scalar W1 = a(m - 1, m - 1) / b(m - 1, m - 1) - a(1, 1) / b(1, 1), W2 = a(m, m) / b(m, m) - a(1, 1) / b(1, 1), + W3 = a(m, m - 1) / b(m - 1, m - 1); + return Vec3( + (W1 * W2 - a(m - 1, m) / b(m, m) * W3 + W3 * b(m - 1, m) / b(m, m) * a(1, 1) / b(1, 1)) * b(1, 1) / a(2, 1) + + a(1, 2) / b(2, 2) - a(1, 1) / b(1, 1) * b(1, 2) / b(2, 2), + (a(2, 2) / b(2, 2) - a(1, 1) / b(1, 1)) - a(2, 1) / b(1, 1) * b(1, 2) / b(2, 2) - W1 - W2 + + W3 * (b(m - 1, m) / b(m, m)), + a(3, 2) / b(2, 2)); +} + +/** \internal Performs one QZ sweep on the unreduced block (S, T)(f:l, f:l). */ +template <typename MatrixType_> +void ComplexQZ<MatrixType_>::do_QZ_step(Index f, Index l, unsigned int iter) { + // Start at the lowest row p > f where the bulge v spreads S(p, p-1) only into negligible entries (the test of + // EISPACK's hqr and LAPACK's xLAHQR, applied to S, with v scaled to |v(0)| + |v(1)| + |v(2)| = 1; rows where that + // sum is 0 or not finite are skipped): + // |S(p, p-1)| (|v(1)| + |v(2)|) <= eps |v(0)| (|S(p-1, p-1)| + |S(p, p)| + |S(p+1, p+1)|). + // From row f the bulge is proportional to S(f+1, f), and a tiny but not negligible entry, e.g. next to zero + // eigenvalues, can stall the iteration. A block of three rows has no such p, as in xLAHQR. + Index p = f; Vec3 X; + for (Index j = l - 2; j > f; --j) { + X = shiftColumn(j, l, iter); + const RealScalar scale = numext::abs(X(0)) + numext::abs(X(1)) + numext::abs(X(2)); + if (!(numext::isfinite)(scale) || scale == RealScalar(0)) continue; + X /= scale; + const RealScalar diagonal = + numext::abs(m_S(j - 1, j - 1)) + numext::abs(m_S(j, j)) + numext::abs(m_S(j + 1, j + 1)); + if (numext::abs(m_S(j, j - 1)) * (numext::abs(X(1)) + numext::abs(X(2))) <= + NumTraits<RealScalar>::epsilon() * numext::abs(X(0)) * diagonal) { + p = j; + break; + } + } + if (p == f) X = shiftColumn(f, l, iter); + const Index m = l - p + 1; // Size of the inner block + Scalar x = X(0), y = X(1), z = X(2); const PermutationMatrix<3, 3, int> S3(Vector3i(2, 0, 1)); for (Index k = p; k < p + m - 2; k++) { X << x, y, z; @@ -517,6 +542,11 @@ m_S.template middleRows<3>(k) .rightCols((std::min)(m_n, m_n - k + 1)) .applyHouseholderOnTheLeft(ess, tau, m_ws.data()); + if (k > f && k == p) { + // The reflector spread the negligible S(p, p-1) into the entries below it, which are dropped. + m_S.coeffRef(k + 1, k - 1) = Scalar(0); + m_S.coeffRef(k + 2, k - 1) = Scalar(0); + } m_T.template middleRows<3>(k).rightCols(m_n - k).applyHouseholderOnTheLeft(ess, tau, m_ws.data()); if (m_computeQZ) m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(ess, numext::conj(tau), m_ws.data()); @@ -612,19 +642,73 @@ if (m_computeQZ) m_Z.applyOnTheLeft(i, i + 1, G); - Mat2 tildeSi = m_S.template block<2, 2>(i, i), tildeTi = m_T.template block<2, 2>(i, i); - Mat2 C = tildeSi.norm() < (lambda * tildeTi).norm() ? tildeSi : lambda * tildeTi; - G.makeGivens(C(0, 0), C(1, 0)); + // The left rotation from T's first column keeps T triangular; the residual ~eps (|S| + |lambda| |T|) of the null + // vector lands in S(i+1, i), where the refinement below removes it. + if (numext::abs(m_T(i, i)) + numext::abs(m_T(i + 1, i)) > tolT) + G.makeGivens(m_T(i, i), m_T(i + 1, i)); + else + G.makeGivens(m_S(i, i), m_S(i + 1, i)); m_S.applyOnTheLeft(i, i + 1, G.adjoint()); m_T.applyOnTheLeft(i, i + 1, G.adjoint()); if (m_computeQZ) m_Q.applyOnTheRight(i, i + 1, G); } - if (!is_negligible(m_S(i + 1, i), m_normOfS * NumTraits<RealScalar>::epsilon())) { + // Rotations formed from lambda can leave S(i+1, i) well above eps (|S(i, i)| + |S(i+1, i+1)|), e.g. for graded T. + // Refine with single-shift QZ steps. With K = T^-1 S = [a b; c d], the Rayleigh quotient shift d = S(i+1, i+1) / + // T(i+1, i+1) maps |c| to |c|^2 |b| / (|a - d|^2 + |c|^2): quadratic for separated eigenvalues, but not monotone. Up + // to four steps run on a copy of the block. Of the prefixes that pass the final test |S(i+1, i)| <= eps m_normOfS, or + // of all if none does, the one with the smallest |S(i+1, i)| / (|S(i, i)| + |S(i+1, i+1)|) is applied. + const auto relativeSubdiagonal = [](const Mat2& s) { + const RealScalar subdiagonal = numext::abs(s(1, 0)); + const RealScalar diagonal = numext::abs(s(0, 0)) + numext::abs(s(1, 1)); + return subdiagonal <= NumTraits<RealScalar>::epsilon() * diagonal ? RealScalar(0) : subdiagonal / diagonal; + }; + JacobiRotation<Scalar> left[4], right[4]; + Mat2 Sstep = m_S.template block<2, 2>(i, i), Tstep = m_T.template block<2, 2>(i, i); + const RealScalar tolS = m_normOfS * NumTraits<RealScalar>::epsilon(); + RealScalar current = relativeSubdiagonal(Sstep), best = current; + bool bestPasses = is_negligible(Sstep(1, 0), tolS); + int bestSteps = 0; + for (int k = 0; k < 4 && current > RealScalar(0); ++k) { + if (is_negligible(Tstep(1, 1), tolT)) { + // An infinite eigenvalue: deflate it with a rotation of the columns, as in the branch above. + left[k] = JacobiRotation<Scalar>(Scalar(1), Scalar(0)); + right[k].makeGivens(Sstep(1, 1), Sstep(1, 0)); + right[k] = right[k].adjoint(); + } else { + const Scalar shift = Sstep(1, 1) / Tstep(1, 1); + left[k].makeGivens(Sstep(0, 0) - shift * Tstep(0, 0), Sstep(1, 0)); + Sstep.applyOnTheLeft(0, 1, left[k].adjoint()); + Tstep.applyOnTheLeft(0, 1, left[k].adjoint()); + right[k] = computeZk2(Tstep.row(1)); + } + Sstep.applyOnTheRight(0, 1, right[k]); + Tstep.applyOnTheRight(0, 1, right[k]); + Tstep(1, 0) = Scalar(0); + current = relativeSubdiagonal(Sstep); + const bool passes = is_negligible(Sstep(1, 0), tolS); + if ((passes && !bestPasses) || (passes == bestPasses && current < best)) { + best = current; + bestPasses = passes; + bestSteps = k + 1; + } + } + for (int k = 0; k < bestSteps; ++k) { + m_S.rightCols(m_n - i).applyOnTheLeft(i, i + 1, left[k].adjoint()); + m_T.rightCols(m_n - i).applyOnTheLeft(i, i + 1, left[k].adjoint()); + if (m_computeQZ) m_Q.applyOnTheRight(i, i + 1, left[k]); + m_S.topRows(i + 2).applyOnTheRight(i, i + 1, right[k]); + m_T.topRows(i + 2).applyOnTheRight(i, i + 1, right[k]); + m_T(i + 1, i) = Scalar(0); + if (m_computeQZ) m_Z.applyOnTheLeft(i, i + 1, right[k].adjoint()); + } + + if (!is_negligible(m_S(i + 1, i), tolS)) { m_info = ComputationInfo::NumericalIssue; } else { m_S(i + 1, i) = Scalar(0); + m_T(i + 1, i) = Scalar(0); } }
diff --git a/test/complex_qz.cpp b/test/complex_qz.cpp index 8018ca5..20bb66a 100644 --- a/test/complex_qz.cpp +++ b/test/complex_qz.cpp
@@ -60,17 +60,18 @@ template <typename MatrixType, typename QZType> void verify_complex_qz_convergence(const MatrixType& a, const MatrixType& b, const QZType& qz) { using RealScalar = typename MatrixType::RealScalar; - // Accumulated rounding in the unitary transformations is O(n * epsilon). + // Accumulated rounding in the unitary transformations is O(n * epsilon). stableNorm, since the squared entries of + // complex_qz_scaled's pencils underflow. const RealScalar tolerance = RealScalar(128 * a.rows()) * NumTraits<RealScalar>::epsilon(); VERIFY_IS_EQUAL(qz.info(), Success); - VERIFY((a - qz.matrixQ() * qz.matrixS() * qz.matrixZ()).norm() <= tolerance * a.norm()); - VERIFY((b - qz.matrixQ() * qz.matrixT() * qz.matrixZ()).norm() <= tolerance * b.norm()); + VERIFY((a - qz.matrixQ() * qz.matrixS() * qz.matrixZ()).stableNorm() <= tolerance * a.stableNorm()); + VERIFY((b - qz.matrixQ() * qz.matrixT() * qz.matrixZ()).stableNorm() <= tolerance * b.stableNorm()); VERIFY(qz.matrixQ().isUnitary(tolerance)); VERIFY(qz.matrixZ().isUnitary(tolerance)); const MatrixType lowerS = qz.matrixS().template triangularView<StrictlyLower>(); const MatrixType lowerT = qz.matrixT().template triangularView<StrictlyLower>(); - VERIFY(lowerS.norm() <= tolerance * a.norm()); - VERIFY(lowerT.norm() <= tolerance * b.norm()); + VERIFY(lowerS.stableNorm() <= tolerance * a.stableNorm()); + VERIFY(lowerT.stableNorm() <= tolerance * b.stableNorm()); } template <typename MatrixType> @@ -137,6 +138,131 @@ } } +// The DARE example "darex #14" of SciPy's solve_discrete_are tests, from Benner, Laub and Mehrmann, "A Collection of +// Benchmark Examples for the Numerical Solution of Algebraic Riccati Equations II: Discrete-Time Case", TU +// Chemnitz-Zwickau SPC 95_23 (1995): the balanced symplectic pencil with the column of R deflated. Eigenvalues 0 (x3, +// one Jordan chain), infinity (x3) and 1 -+ sqrt(5) 1e-8. With AVX2 and FMA, an iterate had S(1, 0) = 1.8e-17 next to +// diagonal entries 2.4e-15 and 3e-3; sweeps starting at row 0 then stalled for 400 iterations. +template <typename MatrixType> +void complex_qz_darex14() { + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename MatrixType::RealScalar; + const Index m = 4, p = 2 * m + 1; + // H = [A 0 B; -Q I 0; 0 0 R], J = [I 0 0; 0 A^* 0; 0 -B^* 0] with A = (1 - 1e-8) e0 e0^T + shift down, + // B = 1e-8 e0, Q = e3 e3^T and R = 1/4. + MatrixType H = MatrixType::Zero(p, p), J = MatrixType::Zero(p, p); + H(0, 0) = J(m, m) = Scalar(RealScalar(1) - RealScalar(1e-8)); + for (Index i = 1; i < m; ++i) H(i, i - 1) = J(m + i - 1, m + i) = Scalar(1); + H(0, 2 * m) = Scalar(RealScalar(1e-8)); + H(m + 3, 3) = Scalar(-1); + H.block(m, m, m, m).setIdentity(); + H(2 * m, 2 * m) = Scalar(RealScalar(0.25)); + J.topLeftCorner(m, m).setIdentity(); + J(2 * m, m) = Scalar(RealScalar(-1e-8)); + // Symplectic balancing (Benner, SIAM J. Sci. Comput. 22(5), 2001) by diag(D, D^-1, 1/2) with + // D = diag(2^20, 2^14, 2^8, 2^2). + Matrix<RealScalar, Dynamic, 1> scale(p); + const int exponents[m] = {20, 14, 8, 2}; + for (Index i = 0; i < m; ++i) { + scale(i) = std::ldexp(RealScalar(1), exponents[i]); + scale(m + i) = std::ldexp(RealScalar(1), -exponents[i]); + } + scale(2 * m) = RealScalar(0.5); + H = scale.template cast<Scalar>().asDiagonal() * H * scale.cwiseInverse().template cast<Scalar>().asDiagonal(); + J = scale.template cast<Scalar>().asDiagonal() * J * scale.cwiseInverse().template cast<Scalar>().asDiagonal(); + const HouseholderQR<MatrixType> qr(H.rightCols(1)); + const MatrixType a = (qr.householderQ().adjoint() * H.leftCols(2 * m)).bottomRows(2 * m); + const MatrixType b = (qr.householderQ().adjoint() * J.leftCols(2 * m)).bottomRows(2 * m); + + ComplexQZ<MatrixType> qz(a, b); + verify_complex_qz_convergence(a, b, qz); +} + +// a = diag(A1, s A2), b = I with A2 = U diag(lambda) U^* normal and |lambda_k| in [1/2, 1]. Relative deflation +// decomposes the small block as if it were alone, so by Bauer-Fike its eigenvalues move by at most its backward error, +// ~128 n eps s. A normwise deflation threshold eps ||a|| would allow errors ~eps ||A1|| >> s. +template <typename MatrixType> +void complex_qz_reducible(typename MatrixType::RealScalar s) { + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename MatrixType::RealScalar; + using Vector = Matrix<Scalar, Dynamic, 1>; + const Index n1 = 8, n2 = 8, n = n1 + n2; + Vector lambda(n2); + for (Index k = 0; k < n2; ++k) + lambda(k) = std::polar(internal::random<RealScalar>(RealScalar(0.5), RealScalar(1)), + internal::random<RealScalar>(-RealScalar(EIGEN_PI), RealScalar(EIGEN_PI))); + const MatrixType U = HouseholderQR<MatrixType>(MatrixType::Random(n2, n2)).householderQ(); + MatrixType a = MatrixType::Zero(n, n); + a.topLeftCorner(n1, n1).setRandom(); + a.bottomRightCorner(n2, n2) = s * U * lambda.asDiagonal() * U.adjoint(); + const MatrixType b = MatrixType::Identity(n, n); + + ComplexQZ<MatrixType> qz(a, b); + verify_complex_qz_convergence(a, b, qz); + const Vector computed = qz.matrixS().diagonal().cwiseQuotient(qz.matrixT().diagonal()); + const RealScalar tolerance = RealScalar(128 * n) * NumTraits<RealScalar>::epsilon() * s; + for (Index k = 0; k < n2; ++k) VERIFY((computed.array() - s * lambda(k)).abs().minCoeff() <= tolerance); +} + +// Graded pencils leave 2x2 blocks that the closed-form reduction in reduceDiagonal2x2block does not deflate on its own. +template <typename MatrixType> +void complex_qz_graded() { + using RealScalar = typename MatrixType::RealScalar; + for (Index dim : {4, 8}) { + for (int trial = 0; trial < 50; ++trial) { + MatrixType a = MatrixType::Random(dim, dim), b = MatrixType::Random(dim, dim); + for (Index i = 0; i < dim; ++i) { + const RealScalar grade = std::pow(RealScalar(10), -RealScalar(8 * i) / RealScalar(dim)); + a.row(i) *= grade; + b.col(i) *= grade; + } + ComplexQZ<MatrixType> qz(a, b); + verify_complex_qz_convergence(a, b, qz); + } + } +} + +// a = s N for the down-shift N, b = I. At s = 1 the 2x2 blocks have the eigenvalue 0, for which a left rotation formed +// from lambda T is the identity and leaves T(i+1, i) in place. At s = sqrt(min) eps the products of the sweep-start +// test underflow to 0, and an unscaled test accepted the first row it tried (0 <= 0) and dropped fill of size s. +template <typename MatrixType> +void complex_qz_scaled() { + using Scalar = typename MatrixType::Scalar; + using RealScalar = typename MatrixType::RealScalar; + const RealScalar tiny = std::sqrt((std::numeric_limits<RealScalar>::min)()) * NumTraits<RealScalar>::epsilon(); + for (RealScalar s : {RealScalar(1), tiny}) { + for (Index dim : {4, 5, 8}) { + MatrixType a = MatrixType::Zero(dim, dim); + a.diagonal(-1).setConstant(Scalar(s)); + const MatrixType b = MatrixType::Identity(dim, dim); + ComplexQZ<MatrixType> qz(a, b); + verify_complex_qz_convergence(a, b, qz); + } + } +} + +// a = V J W^*, b = V W^* with J the 2x2 Jordan block of the eigenvalue 0. On a defective block, refinement steps in +// reduceDiagonal2x2block can lower |S(1, 0)| / (|S(0, 0)| + |S(1, 1)|) by growing the diagonal while |S(1, 0)| grows +// past eps ||S||. +template <typename MatrixType> +void complex_qz_jordan() { + using Scalar = typename MatrixType::Scalar; + MatrixType J = MatrixType::Zero(2, 2); + J(0, 1) = Scalar(1); + for (int trial = 0; trial < 100; ++trial) { + for (bool real : {true, false}) { + MatrixType V = MatrixType::Random(2, 2), W = MatrixType::Random(2, 2); + if (real) { + V = V.real().template cast<Scalar>(); + W = W.real().template cast<Scalar>(); + } + const MatrixType a = V * J * W.adjoint(), b = V * W.adjoint(); + ComplexQZ<MatrixType> qz(a, b); + verify_complex_qz_convergence(a, b, qz); + } + } +} + // Block Householder operations still allocate at larger sizes (issue #3163). template <typename MatrixType> void complex_qz_no_malloc() { @@ -181,6 +307,24 @@ CALL_SUBTEST_8((complex_qz_exceptional_shift<MatrixXcd>())); CALL_SUBTEST_7((complex_qz_small_b<MatrixXcf>())); CALL_SUBTEST_8((complex_qz_small_b<MatrixXcd>())); + CALL_SUBTEST_8((complex_qz_darex14<MatrixXcd>())); + CALL_SUBTEST_10((complex_qz_darex14<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>())); + CALL_SUBTEST_7((complex_qz_reducible<MatrixXcf>(1e-6f))); + CALL_SUBTEST_8((complex_qz_reducible<MatrixXcd>(1e-10))); + CALL_SUBTEST_9((complex_qz_reducible<Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>>(1e-6f))); + CALL_SUBTEST_10((complex_qz_reducible<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>(1e-10))); + CALL_SUBTEST_7((complex_qz_graded<MatrixXcf>())); + CALL_SUBTEST_8((complex_qz_graded<MatrixXcd>())); + CALL_SUBTEST_9((complex_qz_graded<Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>>())); + CALL_SUBTEST_10((complex_qz_graded<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>())); + CALL_SUBTEST_7((complex_qz_scaled<MatrixXcf>())); + CALL_SUBTEST_8((complex_qz_scaled<MatrixXcd>())); + CALL_SUBTEST_9((complex_qz_scaled<Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>>())); + CALL_SUBTEST_10((complex_qz_scaled<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>())); + CALL_SUBTEST_7((complex_qz_jordan<MatrixXcf>())); + CALL_SUBTEST_8((complex_qz_jordan<MatrixXcd>())); + CALL_SUBTEST_9((complex_qz_jordan<Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>>())); + CALL_SUBTEST_10((complex_qz_jordan<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>())); CALL_SUBTEST_9((complex_qz_exceptional_shift<Matrix<std::complex<float>, Dynamic, Dynamic, RowMajor>>())); CALL_SUBTEST_10((complex_qz_exceptional_shift<Matrix<std::complex<double>, Dynamic, Dynamic, RowMajor>>()));