Skip to content

Commit

Permalink
Compute Eigenvalues for matrix greater than 4x4 numerically
Browse files Browse the repository at this point in the history 
- to calculate the roots of characteristic polynomial is too slow
  Example:
  Eigenvalues({{0,1,1,0,1,0,0,0},{1,0,0,1,0,1,0,0},{1,0,0,1,0,0,1,0},{0,1,1,0,0,0,0,1},{1,0,0,0,0,1,1,0},{0,1,0,0,1,0,0,1},{0,0,1,0,1,0,0,1},{0,0,0,1,0,1,1,0}})
  • Loading branch information
@axkr
axkr committed Sep 22, 2023
1 parent 711d264 commit 20a51e9
Show file tree
Hide file tree
Showing 3 changed files with 68 additions and 24 deletions.
60 changes: 48 additions & 12 deletions ...id_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/LinearAlgebra.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@
import static org.matheclipse.core.expression.F.Subtract;
import static org.matheclipse.core.expression.F.Times;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Comparator;
import java.util.List;
import java.util.function.Predicate;
import org.apache.logging.log4j.LogManager;
Expand Down Expand Up @@ -72,6 +74,7 @@
import org.matheclipse.core.eval.util.IndexTableGenerator;
import org.matheclipse.core.expression.ASTRealMatrix;
import org.matheclipse.core.expression.ASTRealVector;
import org.matheclipse.core.expression.ComplexNum;
import org.matheclipse.core.expression.F;
import org.matheclipse.core.expression.ImplementationStatus;
import org.matheclipse.core.expression.S;
Expand Down Expand Up @@ -2055,17 +2058,19 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
} else {
boolean hasNumericArgument = arg1.isNumericArgument(true);// (IAST.CONTAINS_NUMERIC_ARG);
if (!hasNumericArgument) {
ISymbol x = F.Dummy("x");
IExpr m = engine.evaluate(F.CharacteristicPolynomial(arg1, x));
if (m.isPolynomial(x)) {
IExpr eigenValues =
RootsFunctions.roots(m, false, F.List(x), false, true, engine);
if (eigenValues.isList()) {
if (eigenValues.forAll(v -> v.isNumericFunction())) {
IAST sortFunction =
sortValuesIfNumeric((IASTMutable) eigenValues, numberOfEigenvalues);
if (sortFunction.isPresent()) {
return engine.evaluate(sortFunction);
if (dim[0] <= 4 && dim[1] <= 4) {
ISymbol x = F.Dummy("x");
IExpr m = engine.evaluate(F.CharacteristicPolynomial(arg1, x));
if (m.isPolynomial(x)) {
IExpr eigenValues =
RootsFunctions.roots(m, false, F.List(x), false, true, engine);
if (eigenValues.isList()) {
if (eigenValues.forAll(v -> v.isNumericFunction())) {
IAST sortFunction =
sortValuesIfNumeric((IASTMutable) eigenValues, numberOfEigenvalues);
if (sortFunction.isPresent()) {
return engine.evaluate(sortFunction);
}
}
}
}
Expand All @@ -2081,7 +2086,7 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
// switch to numeric calculation
IExpr eigenValues = numericEval(ast, engine);
if (eigenValues.isList()) {
IAST sortFunction = sortValuesIfNumeric((IASTMutable) eigenValues, numberOfEigenvalues);
IAST sortFunction = pickValues((IASTMutable) eigenValues, numberOfEigenvalues);
if (sortFunction.isPresent()) {
return engine.evaluate(sortFunction);
}
Expand Down Expand Up @@ -2119,11 +2124,41 @@ private IAST sortValuesIfNumeric(IAST eigenValuesList, final IExpr numberOfEigen
return eigenValuesList;
}

private static IAST pickValues(IAST eigenValuesList, final IExpr numberOfEigenvalues) {
if (eigenValuesList.forAll(v -> v.isNumericFunction())) {
if (numberOfEigenvalues != null && numberOfEigenvalues.isPresent()) {
int n = numberOfEigenvalues.toIntDefault();
if (n < 0) {
if (n == Integer.MIN_VALUE) {
return F.NIL;
}
return eigenValuesList.copyFrom(eigenValuesList.size() + n, eigenValuesList.size());
}
return eigenValuesList.copyFrom(1, n + 1);
}
return eigenValuesList;
}
return eigenValuesList;
}

@Override
public IExpr matrixEval(FieldMatrix<IExpr> matrix, Predicate<IExpr> zeroChecker) {
return F.NIL;
}

static final class AbsReverseComparator implements Comparator<Complex> {

@Override
public final int compare(final Complex o1, final Complex o2) {
double n1 = o2.norm();
double n2 = o1.norm();
if (F.isFuzzyEquals(n1, n2, Config.DEFAULT_CHOP_DELTA)) {
return ComplexNum.compare(o1, o2);
}
return n1 < n2 ? -1 : 1;
}
}

@Override
public IAST realMatrixEval(RealMatrix matrix, EvalEngine engine) {
// TODO https://github.com/Hipparchus-Math/hipparchus/issues/174
Expand All @@ -2133,6 +2168,7 @@ public IAST realMatrixEval(RealMatrix matrix, EvalEngine engine) {
ComplexEigenDecomposition.DEFAULT_EPSILON_AV_VD_CHECK, //
(c1, c2) -> Double.compare(c2.norm(), c1.norm()));
Complex[] eigenvalues = ced.getEigenvalues();
Arrays.sort(eigenvalues, 0, eigenvalues.length, new AbsReverseComparator());
return F.mapRange(0, eigenvalues.length, (int i) -> {
if (F.isZero(eigenvalues[i].getImaginary())) {
return F.num(eigenvalues[i].getReal());
Expand Down
28 changes: 16 additions & 12 deletions ...id_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ComplexNum.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -338,38 +338,38 @@ public int compareAbsValueToOne() {
return Double.compare(temp, 1.0);
}

public int compareTo(final Complex that) {
if (fComplex.getReal() < that.getReal()) {
public static int compare(final Complex c1, final Complex c2) {
if (c1.getReal() < c2.getReal()) {
return -1;
}
if (fComplex.getReal() > that.getReal()) {
if (c1.getReal() > c2.getReal()) {
return 1;
}
long l1 = Double.doubleToLongBits(fComplex.getReal());
long l2 = Double.doubleToLongBits(that.getReal());
long l1 = Double.doubleToLongBits(c1.getReal());
long l2 = Double.doubleToLongBits(c2.getReal());
if (l1 < l2) {
return -1;
}
if (l1 > l2) {
return 1;
}

if (F.isZero(that.getImaginary())) {
if (!F.isZero(imDoubleValue())) {
if (F.isZero(c2.getImaginary())) {
if (!F.isZero(c1.getImaginary())) {
return 1;
}
} else if (F.isZero(imDoubleValue())) {
} else if (F.isZero(c1.getImaginary())) {
return -1;
}

if (fComplex.getImaginary() < that.getImaginary()) {
if (c1.getImaginary() < c2.getImaginary()) {
return -1;
}
if (fComplex.getImaginary() > that.getImaginary()) {
if (c1.getImaginary() > c2.getImaginary()) {
return 1;
}
l1 = Double.doubleToLongBits(fComplex.getImaginary());
l2 = Double.doubleToLongBits(that.getImaginary());
l1 = Double.doubleToLongBits(c1.getImaginary());
l2 = Double.doubleToLongBits(c2.getImaginary());
if (l1 < l2) {
return -1;
}
Expand All @@ -379,6 +379,10 @@ public int compareTo(final Complex that) {
return 0;
}

public int compareTo(final Complex that) {
return compare(fComplex, that);
}

/**
* Compares this expression with the specified expression for order. Returns a negative integer,
* zero, or a positive integer as this expression is canonical less than, equal to, or greater
Expand Down
4 changes: 4 additions & 0 deletions ...ary/matheclipse-core/src/test/java/org/matheclipse/core/system/LinearAlgebraTestCase.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -497,6 +497,10 @@ public void testEigensystem() {
}

public void testEigenvalues() {
check(
"Eigenvalues({{0,1,1,0,1,0,0,0},{1,0,0,1,0,1,0,0},{1,0,0,1,0,0,1,0},{0,1,1,0,0,0,0,1},{1,0,0,0,0,1,1,0},{0,1,0,0,1,0,0,1},{0,0,1,0,1,0,0,1},{0,0,0,1,0,1,1,0}})", //
"{-3.0,3.0,-1.0,-1.0,-1.0,1.0,1.0,1.0}");

// print message: Eigenvalues: Sequence specification (+n,-n,{+n},{-n},{m,n}) or {m,n,s}
// expected at position 2 in Eigenvalues({{1,0},{0,1}},{1,1,1,1}).
check("Eigenvalues( {{1,0},{0,1}},{1,1,1,1})", //
Expand Down

0 comments on commit 20a51e9

Please sign in to comment.