Skip to content

Commit 3231855

Browse files
axkraxkr
committed
WIP #979 workaround for singular value exception
- Hipparchus-Math/hipparchus#337 - lower the epsilon criteria for final AV=VD check in Eigen decompositions Config.DEFAULT_EPSILON_AV_VD_CHECK = 1.0e-4
1 parent 2faef51 commit 3231855

File tree

3 files changed

+73
-25
lines changed

3 files changed

+73
-25
lines changed

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/basic/Config.java

+5
Original file line numberDiff line numberDiff line change
@@ -290,6 +290,11 @@ public static Cache<IExpr, Object> getExprCache() {
290290
*/
291291
public static double DEFAULT_ROOTS_CHOP_DELTA = 1.0e-5;
292292

293+
/**
294+
* Epsilon criteria for final AV=VD check in Eigen decompositions.
295+
*/
296+
public static double DEFAULT_EPSILON_AV_VD_CHECK = 1.0e-4;
297+
293298
/**
294299
* Tolerance used in special function algorithms ported from
295300
* <a href="https://github.com/paulmasson/math">math.js</a> and in the JavaScript based plot

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/LinearAlgebra.java

+51-21
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,7 @@
4242
import org.hipparchus.linear.ComplexEigenDecomposition;
4343
import org.hipparchus.linear.DecompositionSolver;
4444
import org.hipparchus.linear.DependentVectorsHandler;
45+
import org.hipparchus.linear.EigenDecompositionNonSymmetric;
4546
import org.hipparchus.linear.FieldDecompositionSolver;
4647
import org.hipparchus.linear.FieldLUDecomposition;
4748
import org.hipparchus.linear.FieldMatrix;
@@ -1986,14 +1987,23 @@ public IExpr matrixEval(FieldMatrix<IExpr> matrix, Predicate<IExpr> zeroChecker)
19861987

19871988
@Override
19881989
public IAST realMatrixEval(RealMatrix matrix, EvalEngine engine) {
1989-
ComplexEigenDecomposition ced = new OrderedComplexEigenDecomposition(matrix, //
1990-
ComplexEigenDecomposition.DEFAULT_EIGENVECTORS_EQUALITY, //
1991-
ComplexEigenDecomposition.DEFAULT_EPSILON, //
1992-
ComplexEigenDecomposition.DEFAULT_EPSILON_AV_VD_CHECK, //
1993-
(c1, c2) -> Double.compare(c2.norm(), c1.norm()));
1994-
IAST eigenvalues = Convert.complexValues2List(ced.getEigenvalues());
1990+
try {
1991+
ComplexEigenDecomposition ced = new OrderedComplexEigenDecomposition(matrix, //
1992+
ComplexEigenDecomposition.DEFAULT_EIGENVECTORS_EQUALITY, //
1993+
ComplexEigenDecomposition.DEFAULT_EPSILON, //
1994+
Config.DEFAULT_EPSILON_AV_VD_CHECK, // ComplexEigenDecomposition.DEFAULT_EPSILON_AV_VD_CHECK
1995+
(c1, c2) -> Double.compare(c2.norm(), c1.norm()));
1996+
IAST eigenvalues = Convert.complexValues2List(ced.getEigenvalues());
1997+
IAST eigenvectors = F.mapRange(0, matrix.getColumnDimension(),
1998+
j -> S.Normalize.of(engine, Convert.complexVector2List(ced.getEigenvector(j))));
1999+
return F.List(eigenvalues, eigenvectors);
2000+
} catch (RuntimeException rex) {
2001+
2002+
}
2003+
EigenDecompositionNonSymmetric edns = new EigenDecompositionNonSymmetric(matrix);
2004+
IAST eigenvalues = Convert.complexValues2List(edns.getEigenvalues());
19952005
IAST eigenvectors = F.mapRange(0, matrix.getColumnDimension(),
1996-
j -> S.Normalize.of(engine, Convert.complexVector2List(ced.getEigenvector(j))));
2006+
j -> S.Normalize.of(engine, Convert.complexVector2List(edns.getEigenvector(j))));
19972007
return F.List(eigenvalues, eigenvectors);
19982008
}
19992009

@@ -2197,13 +2207,26 @@ public IExpr matrixEval(FieldMatrix<IExpr> matrix, Predicate<IExpr> zeroChecker)
21972207

21982208
@Override
21992209
public IAST realMatrixEval(RealMatrix matrix, EvalEngine engine) {
2200-
// TODO https://github.com/Hipparchus-Math/hipparchus/issues/174
2201-
ComplexEigenDecomposition ced = new OrderedComplexEigenDecomposition(matrix, //
2202-
ComplexEigenDecomposition.DEFAULT_EIGENVECTORS_EQUALITY, //
2203-
ComplexEigenDecomposition.DEFAULT_EPSILON, //
2204-
ComplexEigenDecomposition.DEFAULT_EPSILON_AV_VD_CHECK, //
2205-
Comparators.COMPLEX_NORM_REVERSE_COMPARATOR);
2206-
Complex[] eigenvalues = ced.getEigenvalues();
2210+
try {
2211+
// https://github.com/Hipparchus-Math/hipparchus/issues/337
2212+
// TODO https://github.com/Hipparchus-Math/hipparchus/issues/174
2213+
ComplexEigenDecomposition ced = new OrderedComplexEigenDecomposition(matrix, //
2214+
ComplexEigenDecomposition.DEFAULT_EIGENVECTORS_EQUALITY, //
2215+
ComplexEigenDecomposition.DEFAULT_EPSILON, //
2216+
Config.DEFAULT_EPSILON_AV_VD_CHECK, // ComplexEigenDecomposition.DEFAULT_EPSILON_AV_VD_CHECK
2217+
Comparators.COMPLEX_NORM_REVERSE_COMPARATOR);
2218+
Complex[] eigenvalues = ced.getEigenvalues();
2219+
return F.mapRange(0, eigenvalues.length, (int i) -> {
2220+
if (F.isZero(eigenvalues[i].getImaginary())) {
2221+
return F.num(eigenvalues[i].getReal());
2222+
}
2223+
return F.complexNum(eigenvalues[i].getReal(), eigenvalues[i].getImaginary());
2224+
});
2225+
} catch (RuntimeException rex) {
2226+
// rex.printStackTrace();
2227+
}
2228+
EigenDecompositionNonSymmetric edns = new EigenDecompositionNonSymmetric(matrix);
2229+
Complex[] eigenvalues = edns.getEigenvalues();
22072230
return F.mapRange(0, eigenvalues.length, (int i) -> {
22082231
if (F.isZero(eigenvalues[i].getImaginary())) {
22092232
return F.num(eigenvalues[i].getReal());
@@ -2473,14 +2496,21 @@ public IExpr matrixEval(FieldMatrix<IExpr> matrix, Predicate<IExpr> zeroChecker)
24732496

24742497
@Override
24752498
public IAST realMatrixEval(RealMatrix matrix, EvalEngine engine) {
2476-
// TODO https://github.com/Hipparchus-Math/hipparchus/issues/174
2477-
ComplexEigenDecomposition ced = new OrderedComplexEigenDecomposition(matrix, //
2478-
ComplexEigenDecomposition.DEFAULT_EIGENVECTORS_EQUALITY, //
2479-
ComplexEigenDecomposition.DEFAULT_EPSILON, //
2480-
ComplexEigenDecomposition.DEFAULT_EPSILON_AV_VD_CHECK, //
2481-
(c1, c2) -> Double.compare(c2.norm(), c1.norm()));
2499+
try {
2500+
// TODO https://github.com/Hipparchus-Math/hipparchus/issues/174
2501+
ComplexEigenDecomposition ced = new OrderedComplexEigenDecomposition(matrix, //
2502+
ComplexEigenDecomposition.DEFAULT_EIGENVECTORS_EQUALITY, //
2503+
ComplexEigenDecomposition.DEFAULT_EPSILON, //
2504+
Config.DEFAULT_EPSILON_AV_VD_CHECK, // ComplexEigenDecomposition.DEFAULT_EPSILON_AV_VD_CHECK
2505+
(c1, c2) -> Double.compare(c2.norm(), c1.norm()));
2506+
return F.mapRange(0, matrix.getColumnDimension(),
2507+
j -> S.Normalize.of(engine, Convert.complexVector2List(ced.getEigenvector(j))));
2508+
} catch (RuntimeException rex) {
2509+
2510+
}
2511+
EigenDecompositionNonSymmetric edns = new EigenDecompositionNonSymmetric(matrix);
24822512
return F.mapRange(0, matrix.getColumnDimension(),
2483-
j -> S.Normalize.of(engine, Convert.complexVector2List(ced.getEigenvector(j))));
2513+
j -> S.Normalize.of(engine, Convert.complexVector2List(edns.getEigenvector(j))));
24842514
}
24852515

24862516
@Override

‎symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LinearAlgebraTestCase.java

+17-4
Original file line numberDiff line numberDiff line change
@@ -516,6 +516,11 @@ public void testDotIssue932() {
516516

517517
@Test
518518
public void testEigensystem() {
519+
check("Chop(Eigensystem({{1,0,0,0,0},{3,1,0,0,0},{6,3,2,0,0},{10,6,3,2,0},{15,10,6,3,2}}))", //
520+
"{{2.00004,1.99998+I*0.0000311927,1.99998+I*(-0.0000311927),1.0,1.0},{{0,0,0,2.18359*10^-6,1.0},{\n" //
521+
+ "0,0,0,0,0},{0,0,0,0,0},{2.87584*10^-9,0.223607,-0.67082,0.67082,-0.223607},{0,0,\n" //
522+
+ "0,0,0}}}");
523+
519524
// example from https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors
520525
check(//
521526
"Chop(Eigensystem({{2.0,0.0,0.0,0.0},{1.0,2.0,0.0,0.0},{0.0,1.0,3.0,0.0},{0.0,0.0,1.0,3.0}}), 10^-7)", //
@@ -534,6 +539,10 @@ public void testEigensystem() {
534539

535540
@Test
536541
public void testEigenvalues() {
542+
check("Eigenvalues({{1, 0, 0}, {-2, 1, 0}, {0, 1, 1}})", //
543+
"{1,1,1}");
544+
check("Eigenvalues({{1,0,0,0,0},{3,1,0,0,0},{6,3,2,0,0},{10,6,3,2,0},{15,10,6,3,2}})", //
545+
"{2.00004,1.99998+I*(-0.0000311927),1.99998+I*0.0000311927,1.0,1.0}");
537546
check("Eigenvalues({{7}},-1)", //
538547
"{7.0}");
539548
check("Eigenvalues({{-1}},1)", //
@@ -635,6 +644,12 @@ public void testEigenvalues() {
635644

636645
@Test
637646
public void testEigenvectors() {
647+
// TODO https://github.com/Hipparchus-Math/hipparchus/issues/337
648+
check("Eigenvectors({{1, 0, 0}, {-2, 1, 0}, {0, 1, 1}})", //
649+
"{{0.0,0.0,1.0},{0.0,-6.66134*10^-16,1.0},{0.0,6.66134*10^-16,-1.0}}");
650+
check("Eigenvectors({{1,0,0,0,0},{3,1,0,0,0},{6,3,2,0,0},{10,6,3,2,0},{15,10,6,3,2}})", //
651+
"{{0.0,0.0,2.62433*10^-12,2.18359*10^-6,1.0},{0.0,0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0,0.0},{2.87584*10^-9,0.223607,-0.67082,0.67082,-0.223607},{0.0,0.0,0.0,0.0,0.0}}");
652+
638653
check("Eigenvectors(SparseArray({{1.0, 2.0, 3}, {4, 5, 6}, {7, 8, 9}}))", //
639654
"{{0.231971,0.525322,0.818673},{0.78583,0.0867513,-0.612328},{-0.408248,0.816497,-0.408248}}");
640655
check("Eigenvectors({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}) // MatrixForm", //
@@ -676,14 +691,12 @@ public void testEigenvectors() {
676691

677692
@Test
678693
public void testEigenvectorsIssue979() {
694+
// https://github.com/Hipparchus-Math/hipparchus/issues/337
679695
// TODO Hipparchus checks submatrices for singularity
680696
check("Eigenvectors({{1,0,0},\n" //
681697
+ "{-2,1,0},\n"//
682698
+ "{0,1,1}})", //
683-
"Eigenvectors(\n"//
684-
+ "{{1,0,0},\n"//
685-
+ " {-2,1,0},\n"//
686-
+ " {0,1,1}})");
699+
"{{0.0,0.0,1.0},{0.0,-6.66134*10^-16,1.0},{0.0,6.66134*10^-16,-1.0}}");
687700
}
688701

689702
@Test

0 commit comments

Comments
 (0)