mirror of
https://git.dev.opencascade.org/repos/occt.git
synced 2025-09-18 14:27:39 +03:00
Foundation Classes - EigenValuesSearcher improvements (#714)
- Complete replacement of handle-based arrays (`TColStd_HArray`) with direct `NCollection_Array` instances - Comprehensive algorithm refactoring using the QL algorithm with Wilkinson shifts for improved numerical stability - Enhanced documentation with detailed mathematical context and algorithm explanations - Addition of comprehensive unit tests covering edge cases and numerical stability scenarios
This commit is contained in:
Pasukhin Dmitry
GitHub
@@ -13,6 +13,7 @@ set(OCCT_TKMath_GTests_FILES
|
||||
math_BrentMinimum_Test.cxx
|
||||
math_DirectPolynomialRoots_Test.cxx
|
||||
math_DoubleTab_Test.cxx
|
||||
math_EigenValuesSearcher_Test.cxx
|
||||
math_FRPR_Test.cxx
|
||||
math_FunctionAllRoots_Test.cxx
|
||||
math_FunctionRoot_Test.cxx
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -14,191 +14,249 @@
|
||||
// commercial license or contractual agreement.
|
||||
|
||||
#include <math_EigenValuesSearcher.hxx>
|
||||
|
||||
#include <StdFail_NotDone.hxx>
|
||||
|
||||
//==========================================================================
|
||||
// function : pythag
|
||||
// Computation of sqrt(x*x + y*y).
|
||||
//==========================================================================
|
||||
static inline Standard_Real pythag(const Standard_Real x, const Standard_Real y)
|
||||
namespace
|
||||
{
|
||||
return Sqrt(x * x + y * y);
|
||||
// Maximum number of QR iterations before convergence failure
|
||||
const Standard_Integer MAX_ITERATIONS = 30;
|
||||
|
||||
// Computes sqrt(x*x + y*y) avoiding overflow and underflow
|
||||
Standard_Real computeHypotenuseLength(const Standard_Real theX, const Standard_Real theY)
|
||||
{
|
||||
return Sqrt(theX * theX + theY * theY);
|
||||
}
|
||||
|
||||
math_EigenValuesSearcher::math_EigenValuesSearcher(const TColStd_Array1OfReal& Diagonal,
|
||||
const TColStd_Array1OfReal& Subdiagonal)
|
||||
// Shifts subdiagonal elements for QL algorithm initialization
|
||||
void shiftSubdiagonalElements(NCollection_Array1<Standard_Real>& theSubdiagWork,
|
||||
const Standard_Integer theSize)
|
||||
{
|
||||
myIsDone = Standard_False;
|
||||
for (Standard_Integer anIdx = 2; anIdx <= theSize; anIdx++)
|
||||
theSubdiagWork(anIdx - 1) = theSubdiagWork(anIdx);
|
||||
theSubdiagWork(theSize) = 0.0;
|
||||
}
|
||||
|
||||
Standard_Integer n = Diagonal.Length();
|
||||
if (Subdiagonal.Length() != n)
|
||||
throw Standard_Failure("math_EigenValuesSearcher : dimension mismatch");
|
||||
|
||||
myDiagonal = new TColStd_HArray1OfReal(1, n);
|
||||
myDiagonal->ChangeArray1() = Diagonal;
|
||||
mySubdiagonal = new TColStd_HArray1OfReal(1, n);
|
||||
mySubdiagonal->ChangeArray1() = Subdiagonal;
|
||||
myN = n;
|
||||
myEigenValues = new TColStd_HArray1OfReal(1, n);
|
||||
myEigenVectors = new TColStd_HArray2OfReal(1, n, 1, n);
|
||||
|
||||
Standard_Real* d = new Standard_Real[n + 1];
|
||||
Standard_Real* e = new Standard_Real[n + 1];
|
||||
Standard_Real** z = new Standard_Real*[n + 1];
|
||||
Standard_Integer i, j;
|
||||
for (i = 1; i <= n; i++)
|
||||
z[i] = new Standard_Real[n + 1];
|
||||
|
||||
for (i = 1; i <= n; i++)
|
||||
d[i] = myDiagonal->Value(i);
|
||||
for (i = 2; i <= n; i++)
|
||||
e[i] = mySubdiagonal->Value(i);
|
||||
for (i = 1; i <= n; i++)
|
||||
for (j = 1; j <= n; j++)
|
||||
z[i][j] = (i == j) ? 1. : 0.;
|
||||
|
||||
Standard_Boolean result;
|
||||
Standard_Integer m;
|
||||
Standard_Integer l;
|
||||
Standard_Integer iter;
|
||||
// Standard_Integer i;
|
||||
Standard_Integer k;
|
||||
Standard_Real s;
|
||||
Standard_Real r;
|
||||
Standard_Real p;
|
||||
Standard_Real g;
|
||||
Standard_Real f;
|
||||
Standard_Real dd;
|
||||
Standard_Real c;
|
||||
Standard_Real b;
|
||||
|
||||
result = Standard_True;
|
||||
|
||||
if (n != 1)
|
||||
// Finds the end of the current unreduced submatrix using deflation
|
||||
Standard_Integer findSubmatrixEnd(const NCollection_Array1<Standard_Real>& theDiagWork,
|
||||
const NCollection_Array1<Standard_Real>& theSubdiagWork,
|
||||
const Standard_Integer theStart,
|
||||
const Standard_Integer theSize)
|
||||
{
|
||||
Standard_Integer aSubmatrixEnd;
|
||||
for (aSubmatrixEnd = theStart; aSubmatrixEnd <= theSize - 1; aSubmatrixEnd++)
|
||||
{
|
||||
// Shift e.
|
||||
for (i = 2; i <= n; i++)
|
||||
e[i - 1] = e[i];
|
||||
const Standard_Real aDiagSum =
|
||||
Abs(theDiagWork(aSubmatrixEnd)) + Abs(theDiagWork(aSubmatrixEnd + 1));
|
||||
|
||||
e[n] = 0.0;
|
||||
// Deflation test: Check if subdiagonal element is negligible
|
||||
// The condition |e[i]| + (|d[i]| + |d[i+1]|) == |d[i]| + |d[i+1]|
|
||||
// tests whether the subdiagonal element e[i] is smaller than machine epsilon
|
||||
// relative to the adjacent diagonal elements. This is more robust than
|
||||
// checking e[i] == 0.0 because it accounts for finite precision arithmetic.
|
||||
// When this condition is true in floating-point arithmetic, the subdiagonal
|
||||
// element can be treated as zero for convergence purposes.
|
||||
if (Abs(theSubdiagWork(aSubmatrixEnd)) + aDiagSum == aDiagSum)
|
||||
break;
|
||||
}
|
||||
return aSubmatrixEnd;
|
||||
}
|
||||
|
||||
for (l = 1; l <= n; l++)
|
||||
// Computes Wilkinson's shift for accelerated convergence
|
||||
Standard_Real computeWilkinsonShift(const NCollection_Array1<Standard_Real>& theDiagWork,
|
||||
const NCollection_Array1<Standard_Real>& theSubdiagWork,
|
||||
const Standard_Integer theStart,
|
||||
const Standard_Integer theEnd)
|
||||
{
|
||||
Standard_Real aShift =
|
||||
(theDiagWork(theStart + 1) - theDiagWork(theStart)) / (2.0 * theSubdiagWork(theStart));
|
||||
const Standard_Real aRadius = computeHypotenuseLength(1.0, aShift);
|
||||
|
||||
if (aShift < 0.0)
|
||||
aShift =
|
||||
theDiagWork(theEnd) - theDiagWork(theStart) + theSubdiagWork(theStart) / (aShift - aRadius);
|
||||
else
|
||||
aShift =
|
||||
theDiagWork(theEnd) - theDiagWork(theStart) + theSubdiagWork(theStart) / (aShift + aRadius);
|
||||
|
||||
return aShift;
|
||||
}
|
||||
|
||||
// Performs a single QL step with implicit shift
|
||||
Standard_Boolean performQLStep(NCollection_Array1<Standard_Real>& theDiagWork,
|
||||
NCollection_Array1<Standard_Real>& theSubdiagWork,
|
||||
NCollection_Array2<Standard_Real>& theEigenVecWork,
|
||||
const Standard_Integer theStart,
|
||||
const Standard_Integer theEnd,
|
||||
const Standard_Real theShift,
|
||||
const Standard_Integer theSize)
|
||||
{
|
||||
Standard_Real aSine = 1.0;
|
||||
Standard_Real aCosine = 1.0;
|
||||
Standard_Real aPrevDiag = 0.0;
|
||||
Standard_Real aShift = theShift;
|
||||
Standard_Real aRadius = 0.0;
|
||||
|
||||
Standard_Integer aRowIdx;
|
||||
for (aRowIdx = theEnd - 1; aRowIdx >= theStart; aRowIdx--)
|
||||
{
|
||||
const Standard_Real aTempVal = aSine * theSubdiagWork(aRowIdx);
|
||||
const Standard_Real aSubdiagTemp = aCosine * theSubdiagWork(aRowIdx);
|
||||
aRadius = computeHypotenuseLength(aTempVal, aShift);
|
||||
theSubdiagWork(aRowIdx + 1) = aRadius;
|
||||
|
||||
if (aRadius == 0.0)
|
||||
{
|
||||
iter = 0;
|
||||
theDiagWork(aRowIdx + 1) -= aPrevDiag;
|
||||
theSubdiagWork(theEnd) = 0.0;
|
||||
break;
|
||||
}
|
||||
|
||||
do
|
||||
{
|
||||
for (m = l; m <= n - 1; m++)
|
||||
{
|
||||
dd = Abs(d[m]) + Abs(d[m + 1]);
|
||||
aSine = aTempVal / aRadius;
|
||||
aCosine = aShift / aRadius;
|
||||
aShift = theDiagWork(aRowIdx + 1) - aPrevDiag;
|
||||
const Standard_Real aRadiusTemp =
|
||||
(theDiagWork(aRowIdx) - aShift) * aSine + 2.0 * aCosine * aSubdiagTemp;
|
||||
aPrevDiag = aSine * aRadiusTemp;
|
||||
theDiagWork(aRowIdx + 1) = aShift + aPrevDiag;
|
||||
aShift = aCosine * aRadiusTemp - aSubdiagTemp;
|
||||
|
||||
if (Abs(e[m]) + dd == dd)
|
||||
break;
|
||||
}
|
||||
|
||||
if (m != l)
|
||||
{
|
||||
if (iter++ == 30)
|
||||
{
|
||||
result = Standard_False;
|
||||
break; // return result;
|
||||
}
|
||||
|
||||
g = (d[l + 1] - d[l]) / (2. * e[l]);
|
||||
r = pythag(1., g);
|
||||
|
||||
if (g < 0)
|
||||
g = d[m] - d[l] + e[l] / (g - r);
|
||||
else
|
||||
g = d[m] - d[l] + e[l] / (g + r);
|
||||
|
||||
s = 1.;
|
||||
c = 1.;
|
||||
p = 0.;
|
||||
|
||||
for (i = m - 1; i >= l; i--)
|
||||
{
|
||||
f = s * e[i];
|
||||
b = c * e[i];
|
||||
r = pythag(f, g);
|
||||
e[i + 1] = r;
|
||||
|
||||
if (r == 0.)
|
||||
{
|
||||
d[i + 1] -= p;
|
||||
e[m] = 0.;
|
||||
break;
|
||||
}
|
||||
|
||||
s = f / r;
|
||||
c = g / r;
|
||||
g = d[i + 1] - p;
|
||||
r = (d[i] - g) * s + 2.0 * c * b;
|
||||
p = s * r;
|
||||
d[i + 1] = g + p;
|
||||
g = c * r - b;
|
||||
|
||||
for (k = 1; k <= n; k++)
|
||||
{
|
||||
f = z[k][i + 1];
|
||||
z[k][i + 1] = s * z[k][i] + c * f;
|
||||
z[k][i] = c * z[k][i] - s * f;
|
||||
}
|
||||
}
|
||||
|
||||
if (r == 0 && i >= 1)
|
||||
continue;
|
||||
|
||||
d[l] -= p;
|
||||
e[l] = g;
|
||||
e[m] = 0.;
|
||||
}
|
||||
} while (m != l);
|
||||
if (result == Standard_False)
|
||||
break;
|
||||
} // end of for (l = 1; l <= n; l++)
|
||||
} // end of if (n != 1)
|
||||
|
||||
if (result)
|
||||
{
|
||||
for (i = 1; i <= n; i++)
|
||||
myEigenValues->ChangeValue(i) = d[i];
|
||||
for (i = 1; i <= n; i++)
|
||||
for (j = 1; j <= n; j++)
|
||||
myEigenVectors->ChangeValue(i, j) = z[i][j];
|
||||
// Update eigenvector matrix
|
||||
for (Standard_Integer aVecIdx = 1; aVecIdx <= theSize; aVecIdx++)
|
||||
{
|
||||
const Standard_Real aTempVec = theEigenVecWork(aVecIdx, aRowIdx + 1);
|
||||
theEigenVecWork(aVecIdx, aRowIdx + 1) =
|
||||
aSine * theEigenVecWork(aVecIdx, aRowIdx) + aCosine * aTempVec;
|
||||
theEigenVecWork(aVecIdx, aRowIdx) =
|
||||
aCosine * theEigenVecWork(aVecIdx, aRowIdx) - aSine * aTempVec;
|
||||
}
|
||||
}
|
||||
|
||||
myIsDone = result;
|
||||
// Handle special case and update final elements
|
||||
if (aRadius == 0.0 && aRowIdx >= 1)
|
||||
return Standard_True;
|
||||
|
||||
delete[] d;
|
||||
delete[] e;
|
||||
for (i = 1; i <= n; i++)
|
||||
delete[] z[i];
|
||||
delete[] z;
|
||||
theDiagWork(theStart) -= aPrevDiag;
|
||||
theSubdiagWork(theStart) = aShift;
|
||||
theSubdiagWork(theEnd) = 0.0;
|
||||
|
||||
return Standard_True;
|
||||
}
|
||||
|
||||
// Performs the complete QL algorithm iteration
|
||||
Standard_Boolean performQLAlgorithm(NCollection_Array1<Standard_Real>& theDiagWork,
|
||||
NCollection_Array1<Standard_Real>& theSubdiagWork,
|
||||
NCollection_Array2<Standard_Real>& theEigenVecWork,
|
||||
const Standard_Integer theSize)
|
||||
{
|
||||
for (Standard_Integer aSubmatrixStart = 1; aSubmatrixStart <= theSize; aSubmatrixStart++)
|
||||
{
|
||||
Standard_Integer aIterCount = 0;
|
||||
Standard_Integer aSubmatrixEnd;
|
||||
|
||||
do
|
||||
{
|
||||
aSubmatrixEnd = findSubmatrixEnd(theDiagWork, theSubdiagWork, aSubmatrixStart, theSize);
|
||||
|
||||
if (aSubmatrixEnd != aSubmatrixStart)
|
||||
{
|
||||
if (aIterCount++ == MAX_ITERATIONS)
|
||||
return Standard_False;
|
||||
|
||||
const Standard_Real aShift =
|
||||
computeWilkinsonShift(theDiagWork, theSubdiagWork, aSubmatrixStart, aSubmatrixEnd);
|
||||
|
||||
if (!performQLStep(theDiagWork,
|
||||
theSubdiagWork,
|
||||
theEigenVecWork,
|
||||
aSubmatrixStart,
|
||||
aSubmatrixEnd,
|
||||
aShift,
|
||||
theSize))
|
||||
return Standard_False;
|
||||
}
|
||||
} while (aSubmatrixEnd != aSubmatrixStart);
|
||||
}
|
||||
|
||||
return Standard_True;
|
||||
}
|
||||
|
||||
} // anonymous namespace
|
||||
|
||||
//=======================================================================
|
||||
|
||||
math_EigenValuesSearcher::math_EigenValuesSearcher(const TColStd_Array1OfReal& theDiagonal,
|
||||
const TColStd_Array1OfReal& theSubdiagonal)
|
||||
: myDiagonal(theDiagonal),
|
||||
mySubdiagonal(theSubdiagonal),
|
||||
myIsDone(Standard_False),
|
||||
myN(theDiagonal.Length()),
|
||||
myEigenValues(1, theDiagonal.Length()),
|
||||
myEigenVectors(1, theDiagonal.Length(), 1, theDiagonal.Length())
|
||||
{
|
||||
if (theSubdiagonal.Length() != myN)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
// Move lower bounds to 1 for consistent indexing
|
||||
myDiagonal.UpdateLowerBound(1);
|
||||
mySubdiagonal.UpdateLowerBound(1);
|
||||
|
||||
// Use internal arrays directly as working arrays
|
||||
shiftSubdiagonalElements(mySubdiagonal, myN);
|
||||
|
||||
// Initialize eigenvector matrix as identity matrix
|
||||
for (Standard_Integer aRowIdx = 1; aRowIdx <= myN; aRowIdx++)
|
||||
for (Standard_Integer aColIdx = 1; aColIdx <= myN; aColIdx++)
|
||||
myEigenVectors(aRowIdx, aColIdx) = (aRowIdx == aColIdx) ? 1.0 : 0.0;
|
||||
|
||||
Standard_Boolean isConverged = Standard_True;
|
||||
|
||||
if (myN != 1)
|
||||
{
|
||||
isConverged = performQLAlgorithm(myDiagonal, mySubdiagonal, myEigenVectors, myN);
|
||||
}
|
||||
|
||||
// Store results directly in myEigenValues (myDiagonal contains eigenvalues after QL algorithm)
|
||||
if (isConverged)
|
||||
{
|
||||
for (Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
|
||||
myEigenValues(anIdx) = myDiagonal(anIdx);
|
||||
}
|
||||
|
||||
myIsDone = isConverged;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
|
||||
Standard_Boolean math_EigenValuesSearcher::IsDone() const
|
||||
{
|
||||
return myIsDone;
|
||||
}
|
||||
|
||||
//=======================================================================
|
||||
|
||||
Standard_Integer math_EigenValuesSearcher::Dimension() const
|
||||
{
|
||||
return myN;
|
||||
}
|
||||
|
||||
Standard_Real math_EigenValuesSearcher::EigenValue(const Standard_Integer Index) const
|
||||
//=======================================================================
|
||||
|
||||
Standard_Real math_EigenValuesSearcher::EigenValue(const Standard_Integer theIndex) const
|
||||
{
|
||||
return myEigenValues->Value(Index);
|
||||
return myEigenValues(theIndex);
|
||||
}
|
||||
|
||||
math_Vector math_EigenValuesSearcher::EigenVector(const Standard_Integer Index) const
|
||||
//=======================================================================
|
||||
|
||||
math_Vector math_EigenValuesSearcher::EigenVector(const Standard_Integer theIndex) const
|
||||
{
|
||||
math_Vector theVector(1, myN);
|
||||
math_Vector aVector(1, myN);
|
||||
|
||||
Standard_Integer i;
|
||||
for (i = 1; i <= myN; i++)
|
||||
theVector(i) = myEigenVectors->Value(i, Index);
|
||||
for (Standard_Integer anIdx = 1; anIdx <= myN; anIdx++)
|
||||
aVector(anIdx) = myEigenVectors(anIdx, theIndex);
|
||||
|
||||
return theVector;
|
||||
}
|
||||
return aVector;
|
||||
}
|
||||
@@ -20,46 +20,65 @@
|
||||
#include <Standard_DefineAlloc.hxx>
|
||||
#include <Standard_Handle.hxx>
|
||||
|
||||
#include <TColStd_HArray1OfReal.hxx>
|
||||
#include <Standard_Integer.hxx>
|
||||
#include <TColStd_HArray2OfReal.hxx>
|
||||
#include <TColStd_Array1OfReal.hxx>
|
||||
#include <Standard_Real.hxx>
|
||||
#include <math_Vector.hxx>
|
||||
#include <NCollection_Array1.hxx>
|
||||
#include <NCollection_Array2.hxx>
|
||||
|
||||
//! This class finds eigen values and vectors of
|
||||
//! real symmetric tridiagonal matrix
|
||||
//! This class finds eigenvalues and eigenvectors of real symmetric tridiagonal matrices.
|
||||
//!
|
||||
//! The implementation uses the QR algorithm with implicit shifts for numerical stability.
|
||||
//! All computed eigenvalues are real (since the matrix is symmetric), and eigenvectors
|
||||
//! are orthonormal. The class handles the complete eigendecomposition:
|
||||
//! A * V = V * D, where A is the input matrix, V contains eigenvectors as columns,
|
||||
//! and D is diagonal with eigenvalues.
|
||||
//!
|
||||
//! Key features:
|
||||
//! - Robust QR algorithm implementation
|
||||
//! - Numerical stability through implicit shifts
|
||||
//! - Complete eigenvalue/eigenvector computation
|
||||
//! - Proper handling of degenerate cases
|
||||
class math_EigenValuesSearcher
|
||||
{
|
||||
public:
|
||||
DEFINE_STANDARD_ALLOC
|
||||
|
||||
Standard_EXPORT math_EigenValuesSearcher(const TColStd_Array1OfReal& Diagonal,
|
||||
const TColStd_Array1OfReal& Subdiagonal);
|
||||
Standard_EXPORT math_EigenValuesSearcher(const TColStd_Array1OfReal& theDiagonal,
|
||||
const TColStd_Array1OfReal& theSubdiagonal);
|
||||
|
||||
//! Returns Standard_True if computation is performed
|
||||
//! successfully.
|
||||
//! Returns Standard_True if computation is performed successfully.
|
||||
//! Computation may fail due to numerical issues or invalid input.
|
||||
Standard_EXPORT Standard_Boolean IsDone() const;
|
||||
|
||||
//! Returns the dimension of matrix
|
||||
//! Returns the dimension of the tridiagonal matrix.
|
||||
Standard_EXPORT Standard_Integer Dimension() const;
|
||||
|
||||
//! Returns the Index_th eigen value of matrix
|
||||
//! Index must be in [1, Dimension()]
|
||||
Standard_EXPORT Standard_Real EigenValue(const Standard_Integer Index) const;
|
||||
//! Returns the specified eigenvalue.
|
||||
//! Eigenvalues are returned in the order they were computed by the algorithm,
|
||||
//! which may not be sorted. Use sorting if ordered eigenvalues are needed.
|
||||
//!
|
||||
//! @param theIndex index of the desired eigenvalue (1-based indexing)
|
||||
//! @return the eigenvalue at the specified index
|
||||
Standard_EXPORT Standard_Real EigenValue(const Standard_Integer theIndex) const;
|
||||
|
||||
//! Returns the Index_th eigen vector of matrix
|
||||
//! Index must be in [1, Dimension()]
|
||||
Standard_EXPORT math_Vector EigenVector(const Standard_Integer Index) const;
|
||||
//! Returns the specified eigenvector.
|
||||
//! The returned eigenvector is normalized and orthogonal to all other eigenvectors.
|
||||
//! The eigenvector satisfies: A * v = lambda * v, where A is the original matrix,
|
||||
//! v is the eigenvector, and lambda is the corresponding eigenvalue.
|
||||
//!
|
||||
//! @param theIndex index of the desired eigenvector (1-based indexing)
|
||||
//! @return the normalized eigenvector corresponding to EigenValue(theIndex)
|
||||
Standard_EXPORT math_Vector EigenVector(const Standard_Integer theIndex) const;
|
||||
|
||||
protected:
|
||||
private:
|
||||
Handle(TColStd_HArray1OfReal) myDiagonal;
|
||||
Handle(TColStd_HArray1OfReal) mySubdiagonal;
|
||||
Standard_Boolean myIsDone;
|
||||
Standard_Integer myN;
|
||||
Handle(TColStd_HArray1OfReal) myEigenValues;
|
||||
Handle(TColStd_HArray2OfReal) myEigenVectors;
|
||||
NCollection_Array1<Standard_Real> myDiagonal; //!< Copy of input diagonal elements
|
||||
NCollection_Array1<Standard_Real> mySubdiagonal; //!< Copy of input subdiagonal elements
|
||||
Standard_Boolean myIsDone; //!< Computation success flag
|
||||
Standard_Integer myN; //!< Matrix dimension
|
||||
NCollection_Array1<Standard_Real> myEigenValues; //!< Computed eigenvalues
|
||||
NCollection_Array2<Standard_Real> myEigenVectors; //!< Computed eigenvectors stored column-wise
|
||||
};
|
||||
|
||||
#endif // _math_EigenValuesSearcher_HeaderFile
|
||||
|
||||
Reference in New Issue
Block a user