org.ojalgo.matrix.decomposition.MatrixDecomposition Maven / Gradle / Ivy
Show all versions of ojalgo Show documentation
/*
* Copyright 1997-2024 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 java.util.function.Supplier;
import org.ojalgo.function.special.MissingMath;
import org.ojalgo.matrix.Provider2D;
import org.ojalgo.matrix.store.ElementsSupplier;
import org.ojalgo.matrix.store.MatrixStore;
import org.ojalgo.matrix.store.PhysicalStore;
import org.ojalgo.matrix.task.DeterminantTask;
import org.ojalgo.matrix.task.InverterTask;
import org.ojalgo.matrix.task.SolverTask;
import org.ojalgo.matrix.transformation.InvertibleFactor;
import org.ojalgo.structure.Access2D;
import org.ojalgo.structure.Access2D.Collectable;
import org.ojalgo.structure.Structure2D;
/**
* Notation used to describe the various matrix decompositions:
*
* - [A] could be any matrix. (The original matrix to decompose.)
* - [A]-1 is the inverse of [A].
* - [A]T is the transpose of [A].
* - [A]H is the conjugate transpose of [A]. [A]H is equilvalent to [A]T if
* the elements are all real.
* - [D] is a diagonal matrix. Possibly bi-, tri- or block-diagonal.
* - [H] is an, upper or lower, Hessenberg matrix.
* - [I] is an identity matrix - obvioulsly orthogonal/unitary.
* - [L] is a lower (left) triangular matrix.
* - [P] is a permutation matrix - an identity matrix with interchanged rows or columns - and
* orthogonal/unitary.
* - [Q] is an orthogonal/unitary matrix. [Q]-1 = [Q]H, and with real matrices = [Q]
* T.
* - [R] is a right (upper) tringular matrix. It is equivalent to [U].
* - [U] is an upper (right) triangular matrix. It is equivalent to [R]. Alternatively [U] is also used to
* denominate the left, orthonormal, singular vectors.
* - [V] is an eigenvector matrix and/or an orthogonal matrix – the columns are the eigenvectors or the
* right, orthonormal, singular vectors.
*
*
* @author apete
*/
public interface MatrixDecomposition> extends Structure2D {
interface Determinant> extends MatrixDecomposition, DeterminantTask, Provider2D.Determinant {
/**
*
* A matrix' determinant is the product of its eigenvalues.
*
*
* @return The matrix' determinant
*/
@Override
N getDeterminant();
@Override
default Provider2D.Determinant toDeterminantProvider(final ElementsSupplier original,
final Supplier> alternativeOriginalSupplier) {
this.decompose(original);
return this;
}
}
/**
* Several matrix decompositions can be expressed "economy sized" - some rows or columns of the decomposed
* matrix parts are not needed for the most releveant use cases, and can therefore be left out. By default
* these matrix decompositions should be "economy sized". Setting {@link #setFullSize(boolean)} to
* true should switch to "full sized".
*
* @author apete
*/
interface EconomySize> extends MatrixDecomposition {
/**
* @return True if it will generate a full sized decomposition.
*/
boolean isFullSize();
}
interface Factory> {
/**
* Will create a new decomposition instance and directly perform the decomposition.
*/
default , DN extends MatrixDecomposition> DN decompose(final Access2D.Collectable> matrix) {
DN retVal = (DN) this.make(matrix);
retVal.decompose(matrix);
return retVal;
}
/**
* Will create a new decomposition instance suitable for typical/most cases.
*
* To calculate the decomposition you then need to call the {@link #decompose(Access2D.Collectable)}
* method.
*
* @return A "decomposer" ready to decompose matrices.
*/
default D make() {
return this.make(TYPICAL);
}
/**
* Will create a new decomposition instance suitable for matrices of the specified size.
*
* To calculate the decomposition you then need to call the {@link #decompose(Access2D.Collectable)}
* method.
*
* @param nbRows The expected number of rows in the matrices to decompose
* @param nbCols The expected number of columns in the matrices to decompose
* @return A "decomposer" ready to decompose matrices.
*/
default D make(final int nbRows, final int nbCols) {
return this.make(new Structure2D() {
@Override
public int getColDim() {
return nbCols;
}
@Override
public int getRowDim() {
return nbRows;
}
});
}
/**
* Will create a new decomposition instance suitable for matrices of the specified size.
*
* To calculate the decomposition you then need to call the {@link #decompose(Access2D.Collectable)}
* method.
*
* @param typical A 2D structure of roughly the expected size with which this decomposition will be
* used.
* @return A "decomposer" ready to decompose matrices.
*/
D make(Structure2D typical);
}
/**
* Some matrix decompositions are only available with hermitian (symmetric) matrices or different
* decomposition algorithms could be used depending on if the matrix is hemitian or not.
*
* @author apete
*/
public interface Hermitian> extends MatrixDecomposition {
/**
* Absolutely must check if the matrix is hermitian or not. Then, depending on the result different
* paths can be chosen - compute or not / choose different algorithms...
*
* @param matrix A matrix to check and then (maybe) decompose
* @return true if the hermitian check passed and decomposition succeeded; false if not
*/
default boolean checkAndDecompose(final MatrixStore matrix) {
this.reset();
if (matrix.isHermitian()) {
return this.decompose(matrix);
} else {
return false;
}
}
}
interface Ordered> extends MatrixDecomposition {
/**
* This is a property of the algorithm/implementation, not the data. Typically relevant for
* {@link SingularValue}, {@link Eigenvalue} or any {@link RankRevealing} decomposition.
*
* @return true if the rows/columns of the returned matrix factors are guaranteed some specific order;
* false if there is no such guarantee.
*/
boolean isOrdered();
}
/**
*
* The pivot or pivot element is the element of a matrix, or an array, which is selected first by an
* algorithm (e.g. Gaussian elimination, simplex algorithm, etc.), to do certain calculations. In the case
* of matrix algorithms, a pivot entry is usually required to be at least distinct from zero, and often
* distant from it; in this case finding this element is called pivoting. Pivoting may be followed by an
* interchange of rows or columns to bring the pivot to a fixed position and allow the algorithm to
* proceed successfully, and possibly to reduce round-off error. It is often used for verifying row
* echelon form.
*
*
* Pivoting might be thought of as swapping or sorting rows or columns in a matrix, and thus it can be
* represented as multiplication by permutation matrices. However, algorithms rarely move the matrix
* elements because this would cost too much time; instead, they just keep track of the permutations.
*
*
* Overall, pivoting adds more operations to the computational cost of an algorithm. These additional
* operations are sometimes necessary for the algorithm to work at all. Other times these additional
* operations are worthwhile because they add numerical stability to the final result.
*
*
* @author apete
*/
interface Pivoting> extends MatrixDecomposition {
/**
* The normal {@link #decompose(Access2D.Collectable)} method must handle cases where pivoting is
* necessary. If you know that pivoting is not needed you may call this method instead - it may be
* faster. Implementing this method, to actually decompose without pivoting, is optional. The default
* implementation simply calls {@link #decompose(Access2D.Collectable)}.
*/
default boolean decomposeWithoutPivoting(final Access2D.Collectable> matrix) {
return this.decompose(matrix);
}
/**
* @return The pivot (row and/or columnn) order
*/
int[] getPivotOrder();
int[] getReversePivotOrder();
/**
* @return true if any pivoting was actually done
*/
boolean isPivoted();
}
/**
* A rank-revealing matrix decomposition of a matrix [A] is a decomposition that is, or can be transformed
* to be, on the form [A]=[X][D][Y]T where:
*
* - [X] and [Y] are square and well conditioned.
* - [D] is diagonal with nonnegative and non-increasing values on the diagonal.
*
*
* The defintion that [X] and [Y] should be well conditioned is subject to interpretation. A specific
* decomposition algorithm can be more or less good at revealing the rank. Typically the
* {@link SingularValue} decomposition is the best.
*
*
* The requirement to have the diagonal elements of [D] ordered can be very practical, but is not always
* strictly necessary in order to just reveal the rank. The method {@link #isOrdered()} indicates if the
* elements (rows and columns) of the returned matrix factors actually are ordered or not for this
* particular implementation.
*
*/
interface RankRevealing> extends Ordered, Provider2D.Rank {
/**
* @param threshold Significance limit
* @return The number of elements in the diagonal matrix that are greater than the threshold
*/
int countSignificant(double threshold);
/**
* The best (and most expensive) way to get the effective numerical rank is by calculating a
* {@link SingularValue} decomposition and then find the number of nonnegligible singular values.
*
* @return The effective numerical rank (best estimate)
*/
@Override
default int getRank() {
return this.countSignificant(this.getRankThreshold());
}
double getRankThreshold();
/**
* @return true if the rank is equal to the minimum of the row and column dimensions; false if not
*/
default boolean isFullRank() {
return this.getRank() == MissingMath.toMinIntExact(this.countRows(), this.countColumns());
}
}
interface Solver> extends MatrixDecomposition, SolverTask, InverterTask, Provider2D.Inverse>>,
Provider2D.Solution>>, InvertibleFactor {
/**
* @param matrix A matrix to decompose
* @return true if the decomposition suceeded AND {@link #isSolvable()}; false if not
*/
default boolean compute(final Collectable> matrix) {
return this.decompose(matrix) && this.isSolvable();
}
@Override
default void ftran(final Collectable> rhs, final PhysicalStore solution) {
this.getSolution(rhs, solution);
}
@Override
default void ftran(final PhysicalStore arg) {
this.ftran(arg, arg);
}
/**
* The output must be a "right inverse" and a "generalised inverse".
*/
MatrixStore getInverse();
/**
*
* Implementing this method is optional.
*
*
* Exactly how a specific implementation makes use of preallocated is not specified by
* this interface. It must be documented for each implementation.
*
*
* Should produce the same results as calling {@link #getInverse()}.
*
*
* @param preallocated Preallocated memory for the results, possibly some intermediate results. You
* must assume this is modified, but you cannot assume it will contain the full/final/correct
* solution. Use {@link #preallocate(int, int)} or {@link #preallocate(Structure2D)} to get a
* suitable instance.
* @return The inverse, this is where you get the solution
* @throws UnsupportedOperationException When/if this feature is not implemented
*/
MatrixStore getInverse(PhysicalStore preallocated);
/**
* [A][X]=[B] or [this][return]=[rhs]
*/
MatrixStore getSolution(Collectable> rhs);
/**
*
* Implementing this method is optional.
*
*
* Exactly how a specific implementation makes use of preallocated is not specified by
* this interface. It must be documented for each implementation.
*
*
* Should produce the same results as calling {@link #getSolution(Collectable)}.
*
*
* @param rhs The Right Hand Side, wont be modfied
* @param preallocated Preallocated memory for the results, possibly some intermediate results. You
* must assume this is modified, but you cannot assume it will contain the full/final/correct
* solution. Use {@link #preallocate(int, int, int)} or
* {@link #preallocate(Structure2D, Structure2D)} to get a suitable instance.
* @return The solution
* @throws UnsupportedOperationException When/if this feature is not implemented
*/
MatrixStore getSolution(Collectable> rhs, PhysicalStore preallocated);
@Override
default Optional> invert() {
if (this.isSolvable()) {
return Optional.of(this.getInverse());
} else {
return Optional.empty();
}
}
/**
* Please note that producing a pseudoinverse and/or a least squares solution is ok! The return value,
* of this method, is not an indication of if the decomposed matrix is square, has full rank, is
* postive definite or whatever. It's that in combination with the specific decomposition algorithm's
* capabilities.
*
* @return true if this matrix decomposition is in a state to be able to deliver an inverse or an
* equation system solution (with some degree of numerical stability).
*/
boolean isSolvable();
@Override
default Optional> solve(final Access2D rhs) {
if (this.isSolvable()) {
return Optional.of(this.getSolution(rhs.asCollectable2D()));
} else {
return Optional.empty();
}
}
@Override
default Provider2D.Inverse>> toInverseProvider(final ElementsSupplier original,
final Supplier> alternativeOriginalSupplier) {
boolean ok = this.decompose(original);
if (ok && this.isSolvable()) {
return this;
} else {
return Optional::empty;
}
}
@Override
default Provider2D.Solution>> toSolutionProvider(final ElementsSupplier body,
final Supplier> alternativeBodySupplier, final Access2D rhs) {
boolean ok = this.decompose(body);
if (ok && this.isSolvable()) {
return this;
} else {
return r -> Optional.empty();
}
}
}
/**
* Eigenvalue and Singular Value decompositions can calculate the "values" only.
*
* @author apete
*/
interface Values> extends Ordered {
/**
* @param matrix The matrix to decompose
* @return The same as {@link Solver#compute(Collectable)} or {@link #decompose(Collectable)} if the
* instance does not implement {@link Solver}.
*/
boolean computeValuesOnly(Access2D.Collectable> matrix);
}
Structure2D TYPICAL = new Structure2D() {
@Override
public int getColDim() {
return 50;
}
@Override
public int getRowDim() {
return 50;
}
};
/**
* @param matrix A matrix to decompose
* @return true if decomposition suceeded; false if not
*/
boolean decompose(Access2D.Collectable> matrix);
/**
* @return true if decomposition has been attemped and was successful; false if not.
* @see #decompose(Access2D.Collectable)
*/
boolean isComputed();
MatrixStore reconstruct();
/**
* Delete computed results, and resets attributes to default values
*/
void reset();
}