diff options
author | Antonio Sanchez <cantonios@google.com> | 2021-03-25 11:08:19 +0000 |
---|---|---|
committer | Antonio Sanchez <cantonios@google.com> | 2021-04-05 11:19:09 -0700 |
commit | 90187a33e141cb8a93953815e8948357350881df (patch) | |
tree | 30029f46bfe5b2be73ad55affd7a9ea9131723e6 /Eigen | |
parent | 3ddc0974ce42d7cdd9161dda2a9558d6800d12c8 (diff) |
Fix SelfAdjoingEigenSolver (#2191)
Adjust the relaxation step to use the condition
```
abs(subdiag[i]) <= epsilon * sqrt(abs(diag[i]) + abs(diag[i+1]))
```
for setting the subdiagonal entry to zero.
Also adjust Wilkinson shift for small `e = subdiag[end-1]` -
I couldn't find a reference for the original, and it was not
consistent with the Wilkinson definition.
Fixes #2191.
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 57 |
1 files changed, 34 insertions, 23 deletions
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 56a12d56f..a6f006c5f 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -498,8 +498,6 @@ template<typename MatrixType, typename DiagType, typename SubDiagType> EIGEN_DEVICE_FUNC ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec) { - EIGEN_USING_STD(abs); - ComputationInfo info; typedef typename MatrixType::Scalar Scalar; @@ -510,15 +508,23 @@ ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag typedef typename DiagType::RealScalar RealScalar; const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); - const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); - + const RealScalar precision_inv = RealScalar(2)/NumTraits<RealScalar>::epsilon(); while (end>0) { - for (Index i = start; i<end; ++i) - if (internal::isMuchSmallerThan(abs(subdiag[i]),(abs(diag[i])+abs(diag[i+1])),precision) || abs(subdiag[i]) <= considerAsZero) - subdiag[i] = 0; + for (Index i = start; i<end; ++i) { + if (numext::abs(subdiag[i]) < considerAsZero) { + subdiag[i] = RealScalar(0); + } else { + // abs(subdiag[i]) <= epsilon * sqrt(abs(diag[i]) + abs(diag[i+1])) + // Scaled to prevent underflows. + const RealScalar scaled_subdiag = precision_inv * subdiag[i]; + if (scaled_subdiag * scaled_subdiag <= (numext::abs(diag[i])+numext::abs(diag[i+1]))) { + subdiag[i] = RealScalar(0); + } + } + } - // find the largest unreduced block + // find the largest unreduced block at the end of the matrix. while (end>0 && subdiag[end-1]==RealScalar(0)) { end--; @@ -821,32 +827,38 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> } namespace internal { + +// Francis implicit QR step. template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> EIGEN_DEVICE_FUNC static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) { - EIGEN_USING_STD(abs); + // Wilkinson Shift. RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); RealScalar e = subdiag[end-1]; // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still // underflow thus leading to inf/NaN values when using the following commented code: -// RealScalar e2 = numext::abs2(subdiag[end-1]); -// RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); + // RealScalar e2 = numext::abs2(subdiag[end-1]); + // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); // This explain the following, somewhat more complicated, version: RealScalar mu = diag[end]; - if(td==RealScalar(0)) - mu -= abs(e); - else - { - RealScalar e2 = numext::abs2(subdiag[end-1]); - RealScalar h = numext::hypot(td,e); - if(e2==RealScalar(0)) mu -= (e / (td + (td>RealScalar(0) ? RealScalar(1) : RealScalar(-1)))) * (e / h); - else mu -= e2 / (td + (td>RealScalar(0) ? h : -h)); + if(td==RealScalar(0)) { + mu -= numext::abs(e); + } else if (e != RealScalar(0)) { + const RealScalar e2 = numext::abs2(e); + const RealScalar h = numext::hypot(td,e); + if(e2 == RealScalar(0)) { + mu -= e / ((td + (td>RealScalar(0) ? h : -h)) / e); + } else { + mu -= e2 / (td + (td>RealScalar(0) ? h : -h)); + } } - + RealScalar x = diag[start] - mu; RealScalar z = subdiag[start]; - for (Index k = start; k < end; ++k) + // If z ever becomes zero, the Givens rotation will be the identity and + // z will stay zero for all future iterations. + for (Index k = start; k < end && z != RealScalar(0); ++k) { JacobiRotation<RealScalar> rot; rot.makeGivens(x, z); @@ -859,12 +871,11 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta diag[k+1] = rot.s() * sdk + rot.c() * dkp1; subdiag[k] = rot.c() * sdk - rot.s() * dkp1; - if (k > start) subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; + // "Chasing the bulge" to return to triangular form. x = subdiag[k]; - if (k < end - 1) { z = -rot.s() * subdiag[k+1]; |