Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

implemented solver as an internal class to avoid building decomposed matrices 

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@728686 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-12-22 13:30:10 +00:00
parent b717e8d1c0
commit bdf20e5cf9
3 changed files with 224 additions and 110 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/java/org/apache/commons/math/linear/EigenDecomposition.java
View File

@ -37,6 +37,8 @@ import java.io.Serializable;
* <li>the <code>getRealEigenvalues</code> method has been renamed as {@link
* #getEigenValues() getEigenValues},</li>
* <li>the <code>getImagEigenvalues</code> method has been removed</li>
* <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
* <li>a {@link #getSolver() getSolver} method has been added.</li>
* </ul>
* @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
* @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
@ -93,4 +95,16 @@ public interface EigenDecomposition extends Serializable {
*/
RealVector getEigenvector(int i);
/**
* Return the determinant of the matrix
* @return determinant of the matrix
*/
double getDeterminant();
/**
* Get a solver for A &times; X = B.
* @return a solver
*/
DecompositionSolver getSolver();
}

198
src/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
View File

@ -55,7 +55,7 @@ import org.apache.commons.math.util.MathUtils;
public class EigenDecompositionImpl implements EigenDecomposition {
/** Serializable version identifier. */
private static final long serialVersionUID = 3125911889630623276L;
private static final long serialVersionUID = 1625101476333719659L;
/** Tolerance. */
private static final double TOLERANCE = 100 * MathUtils.EPSILON;
@ -300,6 +300,202 @@ public class EigenDecompositionImpl implements EigenDecomposition {
return eigenvectors[i].copy();
}
/**
* Return the determinant of the matrix
* @return determinant of the matrix
* @see #isNonSingular()
*/
public double getDeterminant() {
double determinant = 1;
for (double lambda : eigenvalues) {
determinant *= lambda;
}
return determinant;
}
/** {@inheritDoc} */
public DecompositionSolver getSolver() {
if (eigenvectors == null) {
findEigenVectors();
}
return new Solver(eigenvalues, eigenvectors);
}
/** Specialized solver. */
private static class Solver implements DecompositionSolver {
/** Serializable version identifier. */
private static final long serialVersionUID = -8965845906036558410L;
/** Eigenvalues. */
private final double[] eigenvalues;
/** Eigenvectors. */
private final RealVectorImpl[] eigenvectors;
/**
* Build a solver from decomposed matrix.
* @param eigenvalues eigenvalues
* @param eigenvectors eigenvectors
*/
private Solver(final double[] eigenvalues, final RealVectorImpl[] eigenvectors) {
this.eigenvalues = eigenvalues;
this.eigenvectors = eigenvectors;
}
/** Solve the linear equation A &times; X = B for symmetric matrices A.
* <p>This method only find exact linear solutions, i.e. solutions for
* which ||A &times; X - B|| is exactly 0.</p>
* @param b right-hand side of the equation A &times; X = B
* @return a vector X that minimizes the two norm of A &times; X - B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
*/
public double[] solve(final double[] b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = eigenvalues.length;
if (b.length != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] bp = new double[m];
for (int i = 0; i < m; ++i) {
final RealVectorImpl v = eigenvectors[i];
final double[] vData = v.getDataRef();
final double s = v.dotProduct(b) / eigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j] += s * vData[j];
}
}
return bp;
}
/** Solve the linear equation A &times; X = B for symmetric matrices A.
* <p>This method only find exact linear solutions, i.e. solutions for
* which ||A &times; X - B|| is exactly 0.</p>
* @param b right-hand side of the equation A &times; X = B
* @return a vector X that minimizes the two norm of A &times; X - B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
*/
public RealVector solve(final RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = eigenvalues.length;
if (b.getDimension() != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] bp = new double[m];
for (int i = 0; i < m; ++i) {
final RealVectorImpl v = eigenvectors[i];
final double[] vData = v.getDataRef();
final double s = v.dotProduct(b) / eigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j] += s * vData[j];
}
}
return new RealVectorImpl(bp, false);
}
/** Solve the linear equation A &times; X = B for symmetric matrices A.
* <p>This method only find exact linear solutions, i.e. solutions for
* which ||A &times; X - B|| is exactly 0.</p>
* @param b right-hand side of the equation A &times; X = B
* @return a matrix X that minimizes the two norm of A &times; X - B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix solve(final RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = eigenvalues.length;
if (b.getRowDimension() != m) {
throw new IllegalArgumentException("Incorrect row dimension");
}
final int nColB = b.getColumnDimension();
final double[][] bp = new double[m][nColB];
for (int k = 0; k < nColB; ++k) {
for (int i = 0; i < m; ++i) {
final RealVectorImpl v = eigenvectors[i];
final double[] vData = v.getDataRef();
double s = 0;
for (int j = 0; j < m; ++j) {
s += v.getEntry(j) * b.getEntry(j, k);
}
s /= eigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j][k] += s * vData[j];
}
}
}
return MatrixUtils.createRealMatrix(bp);
}
/**
* Check if the decomposed matrix is non-singular.
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular() {
for (double lambda : eigenvalues) {
if (lambda == 0) {
return false;
}
}
return true;
}
/** Get the inverse of the decomposed matrix.
* @return inverse matrix
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix getInverse()
throws InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = eigenvalues.length;
final double[][] invData = new double[m][m];
for (int i = 0; i < m; ++i) {
final double[] invI = invData[i];
for (int j = 0; j < m; ++j) {
double invIJ = 0;
for (int k = 0; k < m; ++k) {
final double[] vK = eigenvectors[k].getDataRef();
invIJ += vK[i] * vK[j] / eigenvalues[k];
}
invI[j] = invIJ;
}
}
return MatrixUtils.createRealMatrix(invData);
}
}
/**
* Transform matrix to tridiagonal.
* @param matrix matrix to transform

122
src/java/org/apache/commons/math/linear/EigenSolver.java
View File

@ -29,17 +29,21 @@ package org.apache.commons.math.linear;
public class EigenSolver implements DecompositionSolver {
/** Serializable version identifier. */
private static final long serialVersionUID = 4339008311386325953L;
private static final long serialVersionUID = -74798755223915020L;
/** Underlying decomposition. */
private final EigenDecomposition decomposition;
/** Underlying solver. */
private final DecompositionSolver solver;
/** Determinant. */
private final double determinant;
/**
* Simple constructor.
* @param decomposition decomposition to use
*/
public EigenSolver(final EigenDecomposition decomposition) {
this.decomposition = decomposition;
this.solver = decomposition.getSolver();
this.determinant = decomposition.getDeterminant();
}
/** Solve the linear equation A &times; X = B for symmetric matrices A.
@ -52,28 +56,7 @@ public class EigenSolver implements DecompositionSolver {
*/
public double[] solve(final double[] b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final double[] eigenvalues = decomposition.getEigenvalues();
final int m = eigenvalues.length;
if (b.length != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] bp = new double[m];
for (int i = 0; i < m; ++i) {
final RealVector v = decomposition.getEigenvector(i);
final double s = v.dotProduct(b) / eigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j] += s * v.getEntry(j);
}
}
return bp;
return solver.solve(b);
}
/** Solve the linear equation A &times; X = B for symmetric matrices A.
@ -86,28 +69,7 @@ public class EigenSolver implements DecompositionSolver {
*/
public RealVector solve(final RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final double[] eigenvalues = decomposition.getEigenvalues();
final int m = eigenvalues.length;
if (b.getDimension() != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] bp = new double[m];
for (int i = 0; i < m; ++i) {
final RealVector v = decomposition.getEigenvector(i);
final double s = v.dotProduct(b) / eigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j] += s * v.getEntry(j);
}
}
return new RealVectorImpl(bp, false);
return solver.solve(b);
}
/** Solve the linear equation A &times; X = B for symmetric matrices A.
@ -120,35 +82,7 @@ public class EigenSolver implements DecompositionSolver {
*/
public RealMatrix solve(final RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final double[] eigenvalues = decomposition.getEigenvalues();
final int m = eigenvalues.length;
if (b.getRowDimension() != m) {
throw new IllegalArgumentException("Incorrect row dimension");
}
final int nColB = b.getColumnDimension();
final double[][] bp = new double[m][nColB];
for (int k = 0; k < nColB; ++k) {
for (int i = 0; i < m; ++i) {
final RealVector v = decomposition.getEigenvector(i);
double s = 0;
for (int j = 0; j < m; ++j) {
s += v.getEntry(j) * b.getEntry(j, k);
}
s /= eigenvalues[i];
for (int j = 0; j < m; ++j) {
bp[j][k] += s * v.getEntry(j);
}
}
}
return MatrixUtils.createRealMatrix(bp);
return solver.solve(b);
}
/**
@ -157,10 +91,6 @@ public class EigenSolver implements DecompositionSolver {
* @see #isNonSingular()
*/
public double getDeterminant() {
double determinant = 1;
for (double lambda : decomposition.getEigenvalues()) {
determinant *= lambda;
}
return determinant;
}
@ -169,12 +99,7 @@ public class EigenSolver implements DecompositionSolver {
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular() {
for (double lambda : decomposition.getEigenvalues()) {
if (lambda == 0) {
return false;
}
}
return true;
return solver.isNonSingular();
}
/** Get the inverse of the decomposed matrix.
@ -183,28 +108,7 @@ public class EigenSolver implements DecompositionSolver {
*/
public RealMatrix getInverse()
throws InvalidMatrixException {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final double[] eigenvalues = decomposition.getEigenvalues();
final int m = eigenvalues.length;
final double[][] invData = new double[m][m];
for (int i = 0; i < m; ++i) {
final double[] invI = invData[i];
for (int j = 0; j < m; ++j) {
double invIJ = 0;
for (int k = 0; k < m; ++k) {
final RealVector vK = decomposition.getEigenvector(k);
invIJ += vK.getEntry(i) * vK.getEntry(j) / eigenvalues[k];
}
invI[j] = invIJ;
}
}
return MatrixUtils.createRealMatrix(invData);
return solver.getInverse();
}
}