devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix flaky matrix_power test

Browse Source 
libeigen/eigen!2325

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-03-22 09:54:32 -07:00
parent 6490b17e6f
commit ac6aedc60a
3 changed files with 69 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files

17
unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
View File

@@ -178,14 +178,17 @@ void MatrixPowerAtomic<MatrixType>::compute2x2(ResultType& res, RealScalar p) co
for (Index i = 1; i < m_A.cols(); ++i) {
res.coeffRef(i, i) = pow(m_A.coeff(i, i), p);
if (m_A.coeff(i - 1, i - 1) == m_A.coeff(i, i))
res.coeffRef(i - 1, i) = p * pow(m_A.coeff(i, i), p - 1);
else if (2 * abs(m_A.coeff(i - 1, i - 1)) < abs(m_A.coeff(i, i)) ||
2 * abs(m_A.coeff(i, i)) < abs(m_A.coeff(i - 1, i - 1)))
res.coeffRef(i - 1, i) =
(res.coeff(i, i) - res.coeff(i - 1, i - 1)) / (m_A.coeff(i, i) - m_A.coeff(i - 1, i - 1));
Scalar a = m_A.coeff(i - 1, i - 1);
Scalar b = m_A.coeff(i, i);
Scalar diff = b - a;
// Use the derivative formula when eigenvalues are nearly equal to avoid
// catastrophic cancellation in the difference quotient.
if (abs(diff) <= RealScalar(2) * (std::numeric_limits<RealScalar>::epsilon)() * (std::max)(abs(a), abs(b)))
res.coeffRef(i - 1, i) = p * pow(b, p - 1);
else if (2 * abs(a) < abs(b) || 2 * abs(b) < abs(a))
res.coeffRef(i - 1, i) = (res.coeff(i, i) - res.coeff(i - 1, i - 1)) / diff;
else
res.coeffRef(i - 1, i) = computeSuperDiag(m_A.coeff(i, i), m_A.coeff(i - 1, i - 1), p);
res.coeffRef(i - 1, i) = computeSuperDiag(b, a, p);
res.coeffRef(i - 1, i) *= m_A.coeff(i - 1, i);
}
}

32
unsupported/test/matrix_functions.h
View File

@@ -16,17 +16,41 @@ struct processTriangularMatrix {
static void run(MatrixType&, MatrixType&, const MatrixType&) {}
};
// For real matrices, make sure none of the eigenvalues are negative.
// For real matrices, ensure all eigenvalues have positive real parts
// (needed for matrix log) and cap the condition number.
template <typename MatrixType>
struct processTriangularMatrix<MatrixType, 0> {
typedef typename MatrixType::Scalar Scalar;
static void run(MatrixType& m, MatrixType& T, const MatrixType& U) {
using std::abs;
const Index size = m.cols();
Scalar maxDiag(0);
for (Index i = 0; i < size; ++i) {
if (i == size - 1 || T.coeff(i + 1, i) == 0)
T.coeffRef(i, i) = std::abs(T.coeff(i, i));
else
if (i == size - 1 || numext::is_exactly_zero(T.coeff(i + 1, i))) {
// 1x1 block (real eigenvalue): make positive.
T.coeffRef(i, i) = abs(T.coeff(i, i));
} else {
// 2x2 block (complex conjugate pair): eigenvalues are T(i,i) ± bi.
// Negate the block if the real part is negative so that the matrix
// log is well-defined (avoids the branch cut on the negative real axis).
if (T.coeff(i, i) < Scalar(0)) {
T.coeffRef(i, i) = -T.coeff(i, i);
T.coeffRef(i + 1, i + 1) = -T.coeff(i + 1, i + 1);
T.coeffRef(i, i + 1) = -T.coeff(i, i + 1);
T.coeffRef(i + 1, i) = -T.coeff(i + 1, i);
}
++i;
}
maxDiag = (std::max)(maxDiag, abs(T.coeff(i, i)));
}
// Clamp small eigenvalues to limit condition number. Matrix power and
// matrix function tests lose too many digits on ill-conditioned matrices.
if (maxDiag > Scalar(0)) {
Scalar minAllowed = maxDiag / Scalar(100);
for (Index i = 0; i < size; ++i) {
if (abs(T.coeff(i, i)) < minAllowed) T.coeffRef(i, i) = minAllowed;
}
}
m = U * T * U.transpose();
}

44
unsupported/test/matrix_power.cpp
View File

@@ -69,6 +69,11 @@ void testGeneral(const MatrixType& m, const typename MatrixType::RealScalar& tol
MatrixType m1, m2, m3, m4, m5;
RealScalar x, y;
// The identities A^(xy) = (A^x)^y and (cA)^y = c^y * A^y involve two
// independent Schur decompositions. Each decomposition + Pade + squaring
// chain introduces rounding that scales with matrix dimension.
RealScalar tol2 = tol * RealScalar(m.rows() * m.rows());
for (int i = 0; i < g_repeat; ++i) {
generateTestMatrix<MatrixType>::run(m1, m.rows());
MatrixPower<MatrixType> mpow(m1);
@@ -84,11 +89,11 @@ void testGeneral(const MatrixType& m, const typename MatrixType::RealScalar& tol
m4 = mpow(x * y);
m5 = m2.pow(y);
VERIFY(m4.isApprox(m5, tol));
VERIFY(m4.isApprox(m5, tol2));
m4 = (std::abs(x) * m1).pow(y);
m5 = std::pow(std::abs(x), y) * m3;
VERIFY(m4.isApprox(m5, tol));
VERIFY(m4.isApprox(m5, tol2));
}
}
@@ -98,11 +103,17 @@ void testSingular(const MatrixType& m_const, const typename MatrixType::RealScal
// MSVC for aligned data types ...
MatrixType& m = const_cast<MatrixType&>(m_const);
typedef typename MatrixType::RealScalar RealScalar;
const int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex;
typedef std::conditional_t<IsComplex, TriangularView<MatrixType, Upper>, const MatrixType&> TriangularType;
std::conditional_t<IsComplex, ComplexSchur<MatrixType>, RealSchur<MatrixType> > schur;
MatrixType T;
// MatrixPower uses its own Schur decomposition while the reference uses
// repeated sqrt of the caller-provided Schur form, so errors compound
// from two independent algorithm paths and scale with matrix dimension.
RealScalar tol2 = tol * RealScalar(m.rows() * m.rows());
for (int i = 0; i < g_repeat; ++i) {
m.setRandom();
m.col(0).fill(0);
@@ -114,13 +125,13 @@ void testSingular(const MatrixType& m_const, const typename MatrixType::RealScal
MatrixPower<MatrixType> mpow(m);
T = T.sqrt();
VERIFY(mpow(0.5L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
VERIFY(mpow(0.5L).isApprox(U * (TriangularType(T) * U.adjoint()), tol2));
T = T.sqrt();
VERIFY(mpow(0.25L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
VERIFY(mpow(0.25L).isApprox(U * (TriangularType(T) * U.adjoint()), tol2));
T = T.sqrt();
VERIFY(mpow(0.125L).isApprox(U * (TriangularType(T) * U.adjoint()), tol));
VERIFY(mpow(0.125L).isApprox(U * (TriangularType(T) * U.adjoint()), tol2));
}
}
@@ -131,12 +142,19 @@ void testLogThenExp(const MatrixType& m_const, const typename MatrixType::RealSc
MatrixType& m = const_cast<MatrixType&>(m_const);
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
Scalar x;
// Computing exp(x*log(m)) chains three separate Pade-based algorithms
// (matrix log, scaling, matrix exp). The accumulated error far exceeds
// what a single matrix power computation achieves.
const Index n = m.rows();
RealScalar tol2 = (std::max)(tol * RealScalar(n * n), test_precision<RealScalar>() * RealScalar(n * n));
for (int i = 0; i < g_repeat; ++i) {
generateTestMatrix<MatrixType>::run(m, m.rows());
x = internal::random<Scalar>();
VERIFY(m.pow(x).isApprox((x * m.log()).exp(), tol));
VERIFY(m.pow(x).isApprox((x * m.log()).exp(), tol2));
}
}
@@ -153,13 +171,13 @@ EIGEN_DECLARE_TEST(matrix_power) {
CALL_SUBTEST_9(test2dHyperbolicRotation<long double>(1e-14L));
CALL_SUBTEST_10(test3dRotation<double>(2e-13));
CALL_SUBTEST_11(test3dRotation<float>(1e-5f));
CALL_SUBTEST_11(test3dRotation<float>(2e-5f));
CALL_SUBTEST_12(test3dRotation<long double>(1e-13L));
CALL_SUBTEST_2(testGeneral(Matrix2d(), 1e-13));
CALL_SUBTEST_2(testGeneral(Matrix2d(), 2e-13));
CALL_SUBTEST_7(testGeneral(Matrix3dRowMajor(), 1e-13));
CALL_SUBTEST_3(testGeneral(Matrix4cd(), 1e-13));
CALL_SUBTEST_4(testGeneral(MatrixXd(8, 8), 2e-12));
CALL_SUBTEST_4(testGeneral(MatrixXd(8, 8), 3e-12));
CALL_SUBTEST_1(testGeneral(Matrix2f(), 1e-4f));
CALL_SUBTEST_5(testGeneral(Matrix3cf(), 1e-4f));
CALL_SUBTEST_8(testGeneral(Matrix4f(), 1e-4f));
@@ -169,10 +187,10 @@ EIGEN_DECLARE_TEST(matrix_power) {
CALL_SUBTEST_11(testGeneral(Matrix3f(), 1e-4f));
CALL_SUBTEST_12(testGeneral(Matrix3e(), 1e-13L));
CALL_SUBTEST_2(testSingular(Matrix2d(), 1e-13));
CALL_SUBTEST_2(testSingular(Matrix2d(), 2e-13));
CALL_SUBTEST_7(testSingular(Matrix3dRowMajor(), 1e-13));
CALL_SUBTEST_3(testSingular(Matrix4cd(), 1e-13));
CALL_SUBTEST_4(testSingular(MatrixXd(8, 8), 2e-12));
CALL_SUBTEST_4(testSingular(MatrixXd(8, 8), 3e-12));
CALL_SUBTEST_1(testSingular(Matrix2f(), 1e-4f));
CALL_SUBTEST_5(testSingular(Matrix3cf(), 1e-4f));
CALL_SUBTEST_8(testSingular(Matrix4f(), 1e-4f));
@@ -182,10 +200,10 @@ EIGEN_DECLARE_TEST(matrix_power) {
CALL_SUBTEST_11(testSingular(Matrix3f(), 1e-4f));
CALL_SUBTEST_12(testSingular(Matrix3e(), 1e-13L));
CALL_SUBTEST_2(testLogThenExp(Matrix2d(), 1e-13));
CALL_SUBTEST_2(testLogThenExp(Matrix2d(), 2e-13));
CALL_SUBTEST_7(testLogThenExp(Matrix3dRowMajor(), 1e-13));
CALL_SUBTEST_3(testLogThenExp(Matrix4cd(), 1e-13));
CALL_SUBTEST_4(testLogThenExp(MatrixXd(8, 8), 2e-12));
CALL_SUBTEST_4(testLogThenExp(MatrixXd(8, 8), 3e-12));
CALL_SUBTEST_1(testLogThenExp(Matrix2f(), 1e-4f));
CALL_SUBTEST_5(testLogThenExp(Matrix3cf(), 1e-4f));
CALL_SUBTEST_8(testLogThenExp(Matrix4f(), 1e-4f));