org.ojalgo.matrix.decomposition.Eigenvalue Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ojalgo Show documentation
Show all versions of ojalgo Show documentation
oj! Algorithms - ojAlgo - is Open Source Java code that has to do with mathematics, linear algebra and optimisation.
/*
* Copyright 1997-2021 Optimatika
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
*/
package org.ojalgo.matrix.decomposition;
import java.util.Optional;
import org.ojalgo.ProgrammingError;
import org.ojalgo.array.Array1D;
import org.ojalgo.array.DenseArray;
import org.ojalgo.matrix.store.GenericStore;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.PhysicalStore;
import org.ojalgo.matrix.store.Primitive64Store;
import org.ojalgo.scalar.ComplexNumber;
import org.ojalgo.scalar.Quaternion;
import org.ojalgo.scalar.RationalNumber;
import org.ojalgo.structure.Access1D;
import org.ojalgo.structure.Access2D;
import org.ojalgo.structure.Structure2D;
import org.ojalgo.type.context.NumberContext;
/**
* [A] = [V][D][V]-1 ([A][V] = [V][D])
*
* - [A] = any square matrix.
* - [V] = contains the eigenvectors as columns.
* - [D] = a diagonal matrix with the eigenvalues on the diagonal (possibly in blocks).
*
*
* [A] is normal if [A][A]H = [A]H[A], and [A] is normal if and only if there exists a
* unitary matrix [Q] such that [A] = [Q][D][Q]H. Hermitian matrices are normal.
*
*
* [V] and [D] can always be calculated in the sense that they will satisfy [A][V] = [V][D], but it is not
* always possible to calculate [V]-1. (Check the rank and/or the condition number of [V] to
* determine the validity of [V][D][V]-1.)
*
*
* The eigenvalues (and their corresponding eigenvectors) of a non-symmetric matrix could be complex.
*
*
* @author apete
*/
public interface Eigenvalue>
extends MatrixDecomposition, MatrixDecomposition.Hermitian, MatrixDecomposition.Determinant, MatrixDecomposition.Values {
public static class Eigenpair implements Comparable {
public final ComplexNumber value;
public final Access1D vector;
public Eigenpair(final ComplexNumber aValue, final Access1D aVector) {
super();
value = aValue;
vector = aVector;
}
public int compareTo(final Eigenpair other) {
return other.value.compareTo(value);
}
@Override
public boolean equals(final Object obj) {
if (this == obj) {
return true;
}
if (obj == null) {
return false;
}
if (!(obj instanceof Eigenpair)) {
return false;
}
final Eigenpair other = (Eigenpair) obj;
if (value == null) {
if (other.value != null) {
return false;
}
} else if (!value.equals(other.value)) {
return false;
}
if (vector == null) {
if (other.vector != null) {
return false;
}
} else if (!vector.equals(other.vector)) {
return false;
}
return true;
}
@Override
public int hashCode() {
final int prime = 31;
int result = 1;
result = (prime * result) + ((value == null) ? 0 : value.hashCode());
result = (prime * result) + ((vector == null) ? 0 : vector.hashCode());
return result;
}
}
interface Factory> extends MatrixDecomposition.Factory> {
default Eigenvalue make(final boolean hermitian) {
return this.make(TYPICAL, hermitian);
}
default Eigenvalue make(final int dimension, final boolean hermitian) {
return this.make(new Structure2D() {
public long countColumns() {
return dimension;
}
public long countRows() {
return dimension;
}
}, hermitian);
}
default Eigenvalue make(final Structure2D typical) {
if (typical instanceof MatrixStore) {
return this.make(typical, ((MatrixStore) typical).isHermitian());
} else {
return this.make(typical, false);
}
}
Eigenvalue make(Structure2D typical, boolean hermitian);
/**
* [A][V] = [B][V][D]
*/
default Eigenvalue.Generalised makeGeneralised(final Structure2D typical) {
return this.makeGeneralised(typical, Eigenvalue.Generalisation.A_B);
}
/**
*
* - http://www.cmth.ph.ic.ac.uk/people/a.mackinnon/Lectures/compphys/node72.html
* - https://www.netlib.org/lapack/lug/node54.html
*
*/
Eigenvalue.Generalised makeGeneralised(Structure2D typical, Eigenvalue.Generalisation type);
}
/**
* Generalized Symmetric Definite
* Eigenproblems
*
* @author apete
*/
enum Generalisation {
/**
* [A][V]=[B][V][D]
*/
A_B,
/**
* [A][B][V]=[V][D]
*/
AB,
/**
* [B][A][V]=[V][D]
*/
BA;
}
interface Generalised> extends Eigenvalue {
/**
* Corresponding to {@link #computeValuesOnly(org.ojalgo.structure.Access2D.Collectable)} but for the
* generalised case.
*
* @see #computeValuesOnly(org.ojalgo.structure.Access2D.Collectable)
*/
default boolean computeValuesOnly(final Access2D.Collectable> matrixA,
final Access2D.Collectable> matrixB) {
return this.prepare(matrixB) && this.computeValuesOnly(matrixA);
}
/**
* Corresponding to {@link #decompose(org.ojalgo.structure.Access2D.Collectable)} but for the
* generalised case.
*
* @see #decompose(org.ojalgo.structure.Access2D.Collectable)
*/
default boolean decompose(final Access2D.Collectable> matrixA,
final Access2D.Collectable> matrixB) {
return this.prepare(matrixB) && this.decompose(matrixA);
}
boolean prepare(Access2D.Collectable> matrixB);
}
Factory COMPLEX = new Factory() {
@Override
public Eigenvalue make(final Structure2D typical, final boolean hermitian) {
return hermitian ? new HermitianEvD.Complex() : null;
}
@Override
public Eigenvalue.Generalised makeGeneralised(final Structure2D typical, final Eigenvalue.Generalisation type) {
PhysicalStore.Factory> factory = GenericStore.COMPLEX;
Cholesky cholesky = Cholesky.COMPLEX.make(typical);
Eigenvalue eigenvalue = this.make(typical, true);
return new GeneralisedEvD<>(factory, cholesky, eigenvalue, type);
}
};
Factory PRIMITIVE = new Factory() {
@Override
public Eigenvalue make(final Structure2D typical) {
if ((8192L < typical.countColumns()) && (typical.count() <= DenseArray.MAX_ARRAY_SIZE)) {
return new DynamicEvD.Primitive();
} else {
return new RawEigenvalue.Dynamic();
}
}
@Override
public Eigenvalue make(final Structure2D typical, final boolean hermitian) {
if (hermitian) {
if ((8192L < typical.countColumns()) && (typical.count() <= DenseArray.MAX_ARRAY_SIZE)) {
return new HermitianEvD.Primitive();
} else {
return new RawEigenvalue.Symmetric();
}
} else {
if ((8192L < typical.countColumns()) && (typical.count() <= DenseArray.MAX_ARRAY_SIZE)) {
return new GeneralEvD.Primitive();
} else {
return new RawEigenvalue.General();
}
}
}
@Override
public Eigenvalue.Generalised makeGeneralised(final Structure2D typical, final Eigenvalue.Generalisation type) {
PhysicalStore.Factory factory = Primitive64Store.FACTORY;
Cholesky cholesky = Cholesky.PRIMITIVE.make(typical);
Eigenvalue eigenvalue = this.make(typical, true);
return new GeneralisedEvD<>(factory, cholesky, eigenvalue, type);
}
};
Factory QUATERNION = new Factory() {
@Override
public Eigenvalue make(final Structure2D typical, final boolean hermitian) {
return hermitian ? new HermitianEvD.Quat() : null;
}
@Override
public Eigenvalue.Generalised makeGeneralised(final Structure2D typical, final Eigenvalue.Generalisation type) {
PhysicalStore.Factory> factory = GenericStore.QUATERNION;
Cholesky cholesky = Cholesky.QUATERNION.make(typical);
Eigenvalue eigenvalue = this.make(typical, true);
return new GeneralisedEvD<>(factory, cholesky, eigenvalue, type);
}
};
Factory RATIONAL = new Factory() {
@Override
public Eigenvalue make(final Structure2D typical, final boolean hermitian) {
return hermitian ? new HermitianEvD.Rational() : null;
}
@Override
public Eigenvalue.Generalised makeGeneralised(final Structure2D typical, final Eigenvalue.Generalisation type) {
PhysicalStore.Factory> factory = GenericStore.RATIONAL;
Cholesky cholesky = Cholesky.RATIONAL.make(typical);
Eigenvalue eigenvalue = this.make(typical, true);
return new GeneralisedEvD<>(factory, cholesky, eigenvalue, type);
}
};
static > boolean equals(final MatrixStore matrix, final Eigenvalue decomposition, final NumberContext context) {
final MatrixStore tmpD = decomposition.getD();
final MatrixStore tmpV = decomposition.getV();
// Check that [A][V] == [V][D] ([A] == [V][D][V]T is not always true)
final MatrixStore tmpStore1 = matrix.multiply(tmpV);
final MatrixStore tmpStore2 = tmpV.multiply(tmpD);
return Access2D.equals(tmpStore1, tmpStore2, context);
}
/**
* @deprecated v48 Use {link #COMPLEX}, {@link #PRIMITIVE}. {@link #QUATERNION} or {@link #RATIONAL}
* innstead.
*/
@Deprecated
@SuppressWarnings("unchecked")
static > Eigenvalue make(final Access2D typical) {
return Eigenvalue.make(typical, Access2D.isHermitian(typical));
}
/**
* @deprecated v48 Use {link #COMPLEX}, {@link #PRIMITIVE}. {@link #QUATERNION} or {@link #RATIONAL}
* innstead.
*/
@Deprecated
@SuppressWarnings("unchecked")
static > Eigenvalue make(final Access2D typical, final boolean hermitian) {
final N tmpNumber = typical.get(0L, 0L);
if (tmpNumber instanceof ComplexNumber) {
return (Eigenvalue) COMPLEX.make(typical, hermitian);
} else if (tmpNumber instanceof Double) {
return (Eigenvalue) PRIMITIVE.make(typical, hermitian);
} else if (tmpNumber instanceof Quaternion) {
return (Eigenvalue) QUATERNION.make(typical, hermitian);
} else if (tmpNumber instanceof RationalNumber) {
return (Eigenvalue) RATIONAL.make(typical, hermitian);
} else {
throw new IllegalArgumentException();
}
}
/**
* @deprecated v48 Use {@link #reconstruct()} instead
*/
@Deprecated
static > MatrixStore reconstruct(final Eigenvalue decomposition) {
return decomposition.reconstruct();
}
/**
* @deprecated With Java 9 this will be made private. Use {@link #getEigenvectors()} or
* {@link #getEigenpair(int)} instead.
*/
@Deprecated
default void copyEigenvector(final int index, final Array1D destination) {
final MatrixStore tmpV = this.getV();
final MatrixStore tmpD = this.getD();
final long tmpDimension = tmpD.countColumns();
final int prevCol = index - 1;
final int nextCol = index + 1;
if ((index < (tmpDimension - 1L)) && (tmpD.doubleValue(nextCol, index) != 0.0)) {
for (int i = 0; i < tmpDimension; i++) {
destination.set(i, ComplexNumber.of(tmpV.doubleValue(i, index), tmpV.doubleValue(i, nextCol)));
}
} else if ((index > 0) && (tmpD.doubleValue(prevCol, index) != 0.0)) {
for (int i = 0; i < tmpDimension; i++) {
destination.set(i, ComplexNumber.of(tmpV.doubleValue(i, prevCol), -tmpV.doubleValue(i, index)));
}
} else {
for (int i = 0; i < tmpDimension; i++) {
destination.set(i, tmpV.doubleValue(i, index));
}
}
}
/**
* The only requirements on [D] are that it should contain the eigenvalues and that [A][V] = [V][D]. The
* ordering of the eigenvalues is not specified.
*
* - If [A] is real and symmetric then [D] is (purely) diagonal with real eigenvalues.
* - If [A] is real but not symmetric then [D] is block-diagonal with real eigenvalues in 1-by-1 blocks
* and complex eigenvalues in 2-by-2 blocks.
* - If [A] is complex then [D] is (purely) diagonal with complex eigenvalues.
*
*
* @return The (block) diagonal eigenvalue matrix.
*/
MatrixStore getD();
default Eigenpair getEigenpair(final int index) {
final long dim = this.getV().countColumns();
final Array1D vector = Array1D.COMPLEX.makeZero(dim);
this.copyEigenvector(index, vector);
final Array1D values = this.getEigenvalues();
final ComplexNumber value = values.get(index);
return new Eigenpair(value, vector);
}
/**
*
* Even for real matrices the eigenvalues (and eigenvectors) are potentially complex numbers. Typically
* they need to be expressed as complex numbers when [A] is not symmetric.
*
*
* The values should be in the same order as the matrices "V" and "D", and if they is ordered or not is
* indicated by the {@link #isOrdered()} method.
*
*
* @return The eigenvalues.
*/
Array1D getEigenvalues();
/**
* @param realParts An array that will receive the real parts of the eigenvalues
* @param imaginaryParts An optional array that, if present, will receive the imaginary parts of the
* eigenvalues
*/
default void getEigenvalues(final double[] realParts, final Optional imaginaryParts) {
ProgrammingError.throwIfNull(realParts, imaginaryParts);
final Array1D values = this.getEigenvalues();
final int length = realParts.length;
if (imaginaryParts.isPresent()) {
final double[] imagParts = imaginaryParts.get();
for (int i = 0; i < length; i++) {
final ComplexNumber value = values.get(i);
realParts[i] = value.getReal();
imagParts[i] = value.getImaginary();
}
} else {
for (int i = 0; i < length; i++) {
realParts[i] = values.doubleValue(i);
}
}
}
/**
* @return A complex valued alternative to {@link #getV()}.
*/
default MatrixStore getEigenvectors() {
final long tmpDimension = this.getV().countColumns();
final GenericStore retVal = GenericStore.COMPLEX.makeZero(tmpDimension, tmpDimension);
for (int j = 0; j < tmpDimension; j++) {
this.copyEigenvector(j, retVal.sliceColumn(0, j));
}
return retVal;
}
// /**
// * @return The matrix exponential
// */
// default MatrixStore getExponential() {
//
// final MatrixStore mtrxV = this.getV();
//
// final PhysicalStore tmpD = this.getD().copy();
// tmpD.modifyDiagonal(mtrxV.physical().function().exp());
// final MatrixStore mtrxD = tmpD.logical().diagonal().get();
//
// return mtrxV.multiply(mtrxD).multiply(mtrxV.conjugate());
// }
//
// /**
// * @return The matrix power
// */
// default MatrixStore getPower(final int exponent) {
//
// final MatrixStore mtrxV = this.getV();
// final MatrixStore mtrxD = this.getD();
//
// MatrixStore retVal = mtrxV;
// for (int e = 0; e < exponent; e++) {
// retVal = retVal.multiply(mtrxD);
// }
// retVal = retVal.multiply(mtrxV.conjugate());
//
// return retVal;
// }
/**
* A matrix' trace is the sum of the diagonal elements. It is also the sum of the eigenvalues. This method
* should return the sum of the eigenvalues.
*
* @return The matrix' trace
*/
ComplexNumber getTrace();
/**
* The columns of [V] represent the eigenvectors of [A] in the sense that [A][V] = [V][D].
*
* @return The eigenvector matrix.
*/
MatrixStore getV();
/**
* If [A] is hermitian then [V][D][V]-1 becomes [Q][D][Q]H...
*/
boolean isHermitian();
/**
* The eigenvalues in D (and the eigenvectors in V) are not necessarily ordered. This is a property of the
* algorithm/implementation, not the data.
*
* @return true if they are ordered
*/
boolean isOrdered();
default MatrixStore reconstruct() {
final MatrixStore mtrxV = this.getV();
MatrixStore mtrxD = this.getD();
return mtrxV.multiply(mtrxD).multiply(mtrxV.conjugate());
}
}