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

improved matrix decomposition API. 

solving a linear system AX = B is now done by a call like:
  RealVector x = new XyzSolver(new XyzDecomposition(a)).solve(b);

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@723736 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-12-05 14:05:50 +00:00
parent cc080433e4
commit 483f55064e
20 changed files with 1593 additions and 1268 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
src/java/org/apache/commons/math/estimation/AbstractEstimator.java
View File

@ -19,8 +19,9 @@ package org.apache.commons.math.estimation;
import java.util.Arrays;
import org.apache.commons.math.linear.DecompositionSolver;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.LUSolver;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealMatrixImpl;
@ -181,8 +182,8 @@ public abstract class AbstractEstimator implements Estimator {
try {
// compute the covariances matrix
DecompositionSolver solver = new DecompositionSolver(new RealMatrixImpl(jTj, false));
RealMatrix inverse = solver.getInverse(solver.luDecompose());
RealMatrix inverse =
new LUSolver(new LUDecompositionImpl(new RealMatrixImpl(jTj, false))).getInverse();
return ((RealMatrixImpl) inverse).getDataRef();
} catch (InvalidMatrixException ime) {
throw new EstimationException("unable to compute covariances: singular problem",

6
src/java/org/apache/commons/math/estimation/GaussNewtonEstimator.java
View File

@ -19,8 +19,9 @@ package org.apache.commons.math.estimation;
import java.io.Serializable;
import org.apache.commons.math.linear.DecompositionSolver;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.LUSolver;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealMatrixImpl;
import org.apache.commons.math.linear.RealVector;
@ -153,8 +154,7 @@ public class GaussNewtonEstimator extends AbstractEstimator implements Serializa
try {
// solve the linearized least squares problem
DecompositionSolver solver = new DecompositionSolver(a);
RealVector dX = solver.solve(b, solver.luDecompose());
RealVector dX = new LUSolver(new LUDecompositionImpl(a)).solve(b);
// update the estimated parameters
for (int i = 0; i < parameters.length; ++i) {

705
src/java/org/apache/commons/math/linear/DecompositionSolver.java
View File

@ -19,15 +19,11 @@ package org.apache.commons.math.linear;
import java.io.Serializable;
import org.apache.commons.math.util.MathUtils;
/**
* Class handling decomposition algorithms that can solve A &times; X = B.
* <p>This class is the entry point for decomposition algorithms like
* {@link QRDecomposition}, {@link LUDecomposition}, {@link
* SingularValueDecomposition} or {@link EigenDecomposition}. All these
* algorithms decompose an A matrix has a product of several specific matrices
* from which they can solve A &times; X = B in least squares sense: they find X
* Interface handling decomposition algorithms that can solve A &times; X = B.
* <p>Decomposition algorithms decompose an A matrix has a product of several specific
* matrices from which they can solve A &times; X = B in least squares sense: they find X
* such that ||A &times; X - B|| is minimal.</p>
* <p>Some solvers like {@link LUDecomposition} can only find the solution for
* square matrices and when the solution is an exact linear solution, i.e. when
@ -38,708 +34,53 @@ import org.apache.commons.math.util.MathUtils;
* @version $Revision$ $Date$
* @since 2.0
*/
public class DecompositionSolver implements Serializable {
public interface DecompositionSolver extends Serializable {
/** Serializable version identifier. */
private static final long serialVersionUID = 182675257956465253L;
/** Matrix to decompose. */
private final RealMatrix matrix;
/**
* Build a decomposition solver for a matrix.
* @param matrix matrix to decompose
*/
public DecompositionSolver(final RealMatrix matrix) {
this.matrix = matrix;
}
/**
* Decompose a matrix using eigendecomposition.
* <p>The split tolerance is set by default to {@link MathUtils#SAFE_MIN}.</p>
* @exception InvalidMatrixException if matrix does not fulfill
* the decomposition requirements (for example non-square matrix
* for {@link LUDecomposition})
*/
public EigenDecomposition eigenDecompose()
throws InvalidMatrixException {
return new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
}
/**
* Decompose a matrix using eigendecomposition.
* @param splitTolerance tolerance on the off-diagonal elements relative to the
* geometric mean to split the tridiagonal matrix (a suggested value is
* {@link MathUtils#SAFE_MIN})
* @exception InvalidMatrixException if matrix does not fulfill
* the decomposition requirements (for example non-square matrix
* for {@link LUDecomposition})
*/
public EigenDecomposition eigenDecompose(final double splitTolerance)
throws InvalidMatrixException {
return new EigenDecompositionImpl(matrix, splitTolerance);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
/** Solve the linear equation A &times; X = B for matrices A.
* <p>The A matrix is implicit, it is provided by the underlying
* decomposition algorithm.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final EigenDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
double[] solve(final double[] b)
throws IllegalArgumentException, InvalidMatrixException;
if (!isNonSingular(decomposition)) {
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;
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
/** Solve the linear equation A &times; X = B for matrices A.
* <p>The A matrix is implicit, it is provided by the underlying
* decomposition algorithm.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final EigenDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
RealVector solve(final RealVector b)
throws IllegalArgumentException, InvalidMatrixException;
if (!isNonSingular(decomposition)) {
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);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
/** Solve the linear equation A &times; X = B for matrices A.
* <p>The A matrix is implicit, it is provided by the underlying
* decomposition algorithm.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final EigenDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
if (!isNonSingular(decomposition)) {
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 new RealMatrixImpl(bp, false);
}
/**
* Return the determinant of the matrix
* @param decomposition decomposition of the matrix A
* @return determinant of the matrix
* @see #isNonSingular()
*/
public double getDeterminant(final EigenDecomposition decomposition) {
double determinant = 1;
for (double lambda : decomposition.getEigenvalues()) {
determinant *= lambda;
}
return determinant;
}
RealMatrix solve(final RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException;
/**
* Check if the decomposed matrix is non-singular.
* @param decomposition decomposition of the matrix A
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular(final EigenDecomposition decomposition) {
for (double lambda : decomposition.getEigenvalues()) {
if (lambda == 0) {
return false;
}
}
return true;
}
boolean isNonSingular();
/** Get the inverse of the decomposed matrix.
/** Get the inverse (or pseudo-inverse) of the decomposed matrix.
* @param decomposition decomposition of the matrix A
* @return inverse matrix
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix getInverse(final EigenDecomposition decomposition)
throws InvalidMatrixException {
if (!isNonSingular(decomposition)) {
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 new RealMatrixImpl(invData, false);
}
/**
* Decompose a matrix using singular value composition.
* @exception InvalidMatrixException if matrix does not fulfill
* the decomposition requirements (for example non-square matrix
* for {@link LUDecomposition})
*/
public SingularValueDecomposition singularDecompose()
throws InvalidMatrixException {
return new SingularValueDecompositionImpl(matrix);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final SingularValueDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
final double[] singularValues = decomposition.getSingularValues();
if (b.length != singularValues.length) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] w = decomposition.getUT().operate(b);
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
w[i] /= si;
}
return decomposition.getV().operate(w);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final SingularValueDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
final double[] singularValues = decomposition.getSingularValues();
if (b.getDimension() != singularValues.length) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final RealVector w = decomposition.getUT().operate(b);
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
w.set(i, w.getEntry(i) / si);
}
return decomposition.getV().operate(w);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final SingularValueDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
final double[] singularValues = decomposition.getSingularValues();
if (b.getRowDimension() != singularValues.length) {
throw new IllegalArgumentException("Incorrect row dimension");
}
final RealMatrixImpl w = (RealMatrixImpl) decomposition.getUT().multiply(b);
final double[][] wData = w.getDataRef();
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
final double inv = 1.0 / si;
final double[] wi = wData[i];
for (int j = 0; j < b.getColumnDimension(); ++j) {
wi[j] *= inv;
}
}
return decomposition.getV().multiply(w);
}
/**
* Check if the decomposed matrix is non-singular.
* @param decomposition decomposition of the matrix A
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular(final SingularValueDecomposition decomposition) {
return decomposition.getRank() == decomposition.getSingularValues().length;
}
/** Get the inverse of the decomposed matrix.
* @param decomposition decomposition of the matrix A
* @return inverse matrix
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix getInverse(final SingularValueDecomposition decomposition)
throws InvalidMatrixException {
if (!isNonSingular(decomposition)) {
throw new SingularMatrixException();
}
return solve(MatrixUtils.createRealIdentityMatrix(decomposition.getSingularValues().length),
decomposition);
}
/**
* Decompose a matrix using QR decomposition.
*/
public QRDecomposition qrDecompose() {
return new QRDecompositionImpl(matrix);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final QRDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
if (decomposition.getR().getRowDimension() != b.length) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (!isNonSingular(decomposition)) {
throw new SingularMatrixException();
}
// solve Q.y = b, using the fact Q is orthogonal
final double[] y = decomposition.getQT().operate(b);
// solve triangular system R.x = y
final RealMatrix r = decomposition.getR();
final double[] x = new double[r.getColumnDimension()];
System.arraycopy(y, 0, x, 0, r.getRowDimension());
for (int i = r.getRowDimension() - 1; i >= 0; --i) {
x[i] /= r.getEntry(i, i);
final double lastX = x[i];
for (int j = i - 1; j >= 0; --j) {
x[j] -= lastX * r.getEntry(j, i);
}
}
return x;
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final QRDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
return new RealVectorImpl(solve(b.getData(), decomposition), false);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final QRDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
if (decomposition.getR().getRowDimension() != b.getRowDimension()) {
throw new IllegalArgumentException("Incorrect row dimension");
}
if (!isNonSingular(decomposition)) {
throw new SingularMatrixException();
}
// solve Q.y = b, using the fact Q is orthogonal
final RealMatrix y = decomposition.getQT().multiply(b);
// solve triangular system R.x = y
final RealMatrix r = decomposition.getR();
final double[][] xData =
new double[r.getColumnDimension()][b.getColumnDimension()];
for (int i = 0; i < r.getRowDimension(); ++i) {
final double[] xi = xData[i];
for (int k = 0; k < xi.length; ++k) {
xi[k] = y.getEntry(i, k);
}
}
for (int i = r.getRowDimension() - 1; i >= 0; --i) {
final double rii = r.getEntry(i, i);
final double[] xi = xData[i];
for (int k = 0; k < xi.length; ++k) {
xi[k] /= rii;
final double lastX = xi[k];
for (int j = i - 1; j >= 0; --j) {
xData[j][k] -= lastX * r.getEntry(j, i);
}
}
}
return new RealMatrixImpl(xData, false);
}
/**
* Check if the decomposed matrix is non-singular.
* @param decomposition decomposition of the matrix A
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular(final QRDecomposition decomposition) {
final RealMatrix r = decomposition.getR();
final int p = Math.min(r.getRowDimension(), r.getColumnDimension());
for (int i = 0; i < p; ++i) {
if (r.getEntry(i, i) == 0) {
return false;
}
}
return true;
}
/** Get the inverse of the decomposed matrix.
* @param decomposition decomposition of the matrix A
* @return inverse matrix
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix getInverse(final QRDecomposition decomposition)
throws InvalidMatrixException {
final RealMatrix r = decomposition.getR();
final int p = Math.min(r.getRowDimension(), r.getColumnDimension());
return solve(MatrixUtils.createRealIdentityMatrix(p), decomposition);
}
/**
* Decompose a matrix using LU decomposition.
* @exception InvalidMatrixException if matrix is non-square)
*/
public LUDecomposition luDecompose()
throws InvalidMatrixException {
return new LUDecompositionImpl(matrix);
}
/**
* Decompose a matrix using LU decomposition.
* @param singularityThreshold threshold (based on partial row norm)
* under which a matrix is considered singular
* @exception InvalidMatrixException if matrix is non-square)
*/
public LUDecomposition luDecompose(final double singularityThreshold)
throws InvalidMatrixException {
return new LUDecompositionImpl(matrix, singularityThreshold);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final LUDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
final int[] pivot = decomposition.getPivot();
final int m = pivot.length;
if (b.length != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (decomposition.isSingular()) {
throw new SingularMatrixException();
}
final double[] bp = new double[m];
// Apply permutations to b
for (int row = 0; row < m; row++) {
bp[row] = b[pivot[row]];
}
// Solve LY = b
final RealMatrix l = decomposition.getL();
for (int col = 0; col < m; col++) {
for (int i = col + 1; i < m; i++) {
bp[i] -= bp[col] * l.getEntry(i, col);
}
}
// Solve UX = Y
final RealMatrix u = decomposition.getU();
for (int col = m - 1; col >= 0; col--) {
bp[col] /= u.getEntry(col, col);
for (int i = 0; i < col; i++) {
bp[i] -= bp[col] * u.getEntry(i, col);
}
}
return bp;
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final LUDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
final int[] pivot = decomposition.getPivot();
final int m = pivot.length;
if (b.getDimension() != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (decomposition.isSingular()) {
throw new SingularMatrixException();
}
final double[] bp = new double[m];
// Apply permutations to b
for (int row = 0; row < m; row++) {
bp[row] = b.getEntry(pivot[row]);
}
// Solve LY = b
final RealMatrix l = decomposition.getL();
for (int col = 0; col < m; col++) {
for (int i = col + 1; i < m; i++) {
bp[i] -= bp[col] * l.getEntry(i, col);
}
}
// Solve UX = Y
final RealMatrix u = decomposition.getU();
for (int col = m - 1; col >= 0; col--) {
bp[col] /= u.getEntry(col, col);
for (int i = 0; i < col; i++) {
bp[i] -= bp[col] * u.getEntry(i, col);
}
}
return new RealVectorImpl(bp, false);
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It <strong>must</strong> have
* already been provided by a previous call to {@link #decompose(RealMatrix)}.</p>
* @param b right-hand side of the equation A &times; X = B
* @param decomposition decomposition of the matrix A
* @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, final LUDecomposition decomposition)
throws IllegalArgumentException, InvalidMatrixException {
final int[] pivot = decomposition.getPivot();
final int m = pivot.length;
if (b.getRowDimension() != m) {
throw new IllegalArgumentException("Incorrect row dimension");
}
if (decomposition.isSingular()) {
throw new SingularMatrixException();
}
final int nColB = b.getColumnDimension();
// Apply permutations to b
final double[][] bp = new double[m][nColB];
for (int row = 0; row < m; row++) {
final double[] bpRow = bp[row];
final int pRow = pivot[row];
for (int col = 0; col < nColB; col++) {
bpRow[col] = b.getEntry(pRow, col);
}
}
// Solve LY = b
final RealMatrix l = decomposition.getL();
for (int col = 0; col < m; col++) {
final double[] bpCol = bp[col];
for (int i = col + 1; i < m; i++) {
final double[] bpI = bp[i];
final double luICol = l.getEntry(i, col);
for (int j = 0; j < nColB; j++) {
bpI[j] -= bpCol[j] * luICol;
}
}
}
// Solve UX = Y
final RealMatrix u = decomposition.getU();
for (int col = m - 1; col >= 0; col--) {
final double[] bpCol = bp[col];
final double luDiag = u.getEntry(col, col);
for (int j = 0; j < nColB; j++) {
bpCol[j] /= luDiag;
}
for (int i = 0; i < col; i++) {
final double[] bpI = bp[i];
final double luICol = u.getEntry(i, col);
for (int j = 0; j < nColB; j++) {
bpI[j] -= bpCol[j] * luICol;
}
}
}
return new RealMatrixImpl(bp, false);
}
/**
* Return the determinant of the matrix
* @param decomposition decomposition of the matrix A
* @return determinant of the matrix
* @see #isNonSingular()
*/
public double getDeterminant(final LUDecomposition decomposition) {
if (decomposition.isSingular()) {
return 0;
} else {
final int m = decomposition.getPivot().length;
final RealMatrix u = decomposition.getU();
double determinant = decomposition.evenPermutation() ? 1 : -1;
for (int i = 0; i < m; i++) {
determinant *= u.getEntry(i, i);
}
return determinant;
}
}
/**
* Check if the decomposed matrix is non-singular.
* @param decomposition decomposition of the matrix A
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular(final LUDecomposition decomposition) {
return !decomposition.isSingular();
}
/** Get the inverse of the decomposed matrix.
* @param decomposition decomposition of the matrix A
* @return inverse matrix
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix getInverse(final LUDecomposition decomposition)
throws InvalidMatrixException {
final int m = decomposition.getPivot().length;
return solve(MatrixUtils.createRealIdentityMatrix(m), decomposition);
}
RealMatrix getInverse()
throws InvalidMatrixException;
}

210
src/java/org/apache/commons/math/linear/EigenSolver.java Normal file
View File

@ -0,0 +1,210 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
/**
* Solver using eigen decomposition to solve A &times; X = B for symmetric matrices A.
* <p>This class finds only exact linear solution, i.e. when
* ||A &times; X - B|| is exactly 0.</p>
*
* @version $Revision$ $Date$
* @since 2.0
*/
public class EigenSolver implements DecompositionSolver {
/** Serializable version identifier. */
private static final long serialVersionUID = 4339008311386325953L;
/** Underlying decomposition. */
private final EigenDecomposition decomposition;
/**
* Simple constructor.
* @param decomposition decomposition to use
*/
public EigenSolver(final EigenDecomposition decomposition) {
this.decomposition = decomposition;
}
/** 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 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;
}
/** 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 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);
}
/** 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 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 new RealMatrixImpl(bp, false);
}
/**
* Return the determinant of the matrix
* @return determinant of the matrix
* @see #isNonSingular()
*/
public double getDeterminant() {
double determinant = 1;
for (double lambda : decomposition.getEigenvalues()) {
determinant *= lambda;
}
return determinant;
}
/**
* Check if the decomposed matrix is non-singular.
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular() {
for (double lambda : decomposition.getEigenvalues()) {
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 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 new RealMatrixImpl(invData, false);
}
}

245
src/java/org/apache/commons/math/linear/LUSolver.java Normal file
View File

@ -0,0 +1,245 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
/**
* Solver using LU decomposition to solve A &times; X = B for square matrices A.
* <p>This class finds only exact linear solution, i.e. when
* ||A &times; X - B|| is exactly 0.</p>
*
* @version $Revision$ $Date$
* @since 2.0
*/
public class LUSolver implements DecompositionSolver {
/** Serializable version identifier. */
private static final long serialVersionUID = -8775006035077527661L;
/** Underlying decomposition. */
private final LUDecomposition decomposition;
/**
* Simple constructor.
* @param decomposition decomposition to use
*/
public LUSolver(final LUDecomposition decomposition) {
this.decomposition = decomposition;
}
/** Solve the linear equation A &times; X = B for square 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 {
final int[] pivot = decomposition.getPivot();
final int m = pivot.length;
if (b.length != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (decomposition.isSingular()) {
throw new SingularMatrixException();
}
final double[] bp = new double[m];
// Apply permutations to b
for (int row = 0; row < m; row++) {
bp[row] = b[pivot[row]];
}
// Solve LY = b
final RealMatrix l = decomposition.getL();
for (int col = 0; col < m; col++) {
for (int i = col + 1; i < m; i++) {
bp[i] -= bp[col] * l.getEntry(i, col);
}
}
// Solve UX = Y
final RealMatrix u = decomposition.getU();
for (int col = m - 1; col >= 0; col--) {
bp[col] /= u.getEntry(col, col);
for (int i = 0; i < col; i++) {
bp[i] -= bp[col] * u.getEntry(i, col);
}
}
return bp;
}
/** Solve the linear equation A &times; X = B for square 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 {
final int[] pivot = decomposition.getPivot();
final int m = pivot.length;
if (b.getDimension() != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (decomposition.isSingular()) {
throw new SingularMatrixException();
}
final double[] bp = new double[m];
// Apply permutations to b
for (int row = 0; row < m; row++) {
bp[row] = b.getEntry(pivot[row]);
}
// Solve LY = b
final RealMatrix l = decomposition.getL();
for (int col = 0; col < m; col++) {
for (int i = col + 1; i < m; i++) {
bp[i] -= bp[col] * l.getEntry(i, col);
}
}
// Solve UX = Y
final RealMatrix u = decomposition.getU();
for (int col = m - 1; col >= 0; col--) {
bp[col] /= u.getEntry(col, col);
for (int i = 0; i < col; i++) {
bp[i] -= bp[col] * u.getEntry(i, col);
}
}
return new RealVectorImpl(bp, false);
}
/** Solve the linear equation A &times; X = B for square 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 {
final int[] pivot = decomposition.getPivot();
final int m = pivot.length;
if (b.getRowDimension() != m) {
throw new IllegalArgumentException("Incorrect row dimension");
}
if (decomposition.isSingular()) {
throw new SingularMatrixException();
}
final int nColB = b.getColumnDimension();
// Apply permutations to b
final double[][] bp = new double[m][nColB];
for (int row = 0; row < m; row++) {
final double[] bpRow = bp[row];
final int pRow = pivot[row];
for (int col = 0; col < nColB; col++) {
bpRow[col] = b.getEntry(pRow, col);
}
}
// Solve LY = b
final RealMatrix l = decomposition.getL();
for (int col = 0; col < m; col++) {
final double[] bpCol = bp[col];
for (int i = col + 1; i < m; i++) {
final double[] bpI = bp[i];
final double luICol = l.getEntry(i, col);
for (int j = 0; j < nColB; j++) {
bpI[j] -= bpCol[j] * luICol;
}
}
}
// Solve UX = Y
final RealMatrix u = decomposition.getU();
for (int col = m - 1; col >= 0; col--) {
final double[] bpCol = bp[col];
final double luDiag = u.getEntry(col, col);
for (int j = 0; j < nColB; j++) {
bpCol[j] /= luDiag;
}
for (int i = 0; i < col; i++) {
final double[] bpI = bp[i];
final double luICol = u.getEntry(i, col);
for (int j = 0; j < nColB; j++) {
bpI[j] -= bpCol[j] * luICol;
}
}
}
return new RealMatrixImpl(bp, false);
}
/**
* Return the determinant of the matrix
* @return determinant of the matrix
* @see #isNonSingular()
*/
public double getDeterminant() {
if (decomposition.isSingular()) {
return 0;
} else {
final int m = decomposition.getPivot().length;
final RealMatrix u = decomposition.getU();
double determinant = decomposition.evenPermutation() ? 1 : -1;
for (int i = 0; i < m; i++) {
determinant *= u.getEntry(i, i);
}
return determinant;
}
}
/**
* Check if the decomposed matrix is non-singular.
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular() {
return !decomposition.isSingular();
}
/** Get the inverse of the decomposed matrix.
* @return inverse matrix
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix getInverse()
throws InvalidMatrixException {
final int m = decomposition.getPivot().length;
return solve(MatrixUtils.createRealIdentityMatrix(m));
}
}

169
src/java/org/apache/commons/math/linear/QRSolver.java Normal file
View File

@ -0,0 +1,169 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
/**
* Class using QR decomposition to solve A &times; X = B in least square sense
* for any matrices A.
* <p>This class solve A &times; X = B in least squares sense: it finds X
* such that ||A &times; X - B|| is minimal.</p>
*
* @version $Revision$ $Date$
* @since 2.0
*/
public class QRSolver implements DecompositionSolver {
/** Serializable version identifier. */
private static final long serialVersionUID = -579465076068393818L;
/** Underlying decomposition. */
private final QRDecomposition decomposition;
/**
* Simple constructor.
* @param decomposition decomposition to use
*/
public QRSolver(final QRDecomposition decomposition) {
this.decomposition = decomposition;
}
/** Solve the linear equation A &times; X = B in least square sense.
* <p>The m&times;n matrix A may not be square, the solution X is
* such that ||A &times; X - B|| is minimal.</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 (decomposition.getR().getRowDimension() != b.length) {
throw new IllegalArgumentException("constant vector has wrong length");
}
if (!isNonSingular()) {
throw new SingularMatrixException();
}
// solve Q.y = b, using the fact Q is orthogonal
final double[] y = decomposition.getQT().operate(b);
// solve triangular system R.x = y
final RealMatrix r = decomposition.getR();
final double[] x = new double[r.getColumnDimension()];
System.arraycopy(y, 0, x, 0, r.getRowDimension());
for (int i = r.getRowDimension() - 1; i >= 0; --i) {
x[i] /= r.getEntry(i, i);
final double lastX = x[i];
for (int j = i - 1; j >= 0; --j) {
x[j] -= lastX * r.getEntry(j, i);
}
}
return x;
}
/** Solve the linear equation A &times; X = B in least square sense.
* <p>The m&times;n matrix A may not be square, the solution X is
* such that ||A &times; X - B|| is minimal.</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 {
return new RealVectorImpl(solve(b.getData()), false);
}
/** Solve the linear equation A &times; X = B in least square sense.
* <p>The m&times;n matrix A may not be square, the solution X is
* such that ||A &times; X - B|| is minimal.</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 (decomposition.getR().getRowDimension() != b.getRowDimension()) {
throw new IllegalArgumentException("Incorrect row dimension");
}
if (!isNonSingular()) {
throw new SingularMatrixException();
}
// solve Q.y = b, using the fact Q is orthogonal
final RealMatrix y = decomposition.getQT().multiply(b);
// solve triangular system R.x = y
final RealMatrix r = decomposition.getR();
final double[][] xData =
new double[r.getColumnDimension()][b.getColumnDimension()];
for (int i = 0; i < r.getRowDimension(); ++i) {
final double[] xi = xData[i];
for (int k = 0; k < xi.length; ++k) {
xi[k] = y.getEntry(i, k);
}
}
for (int i = r.getRowDimension() - 1; i >= 0; --i) {
final double rii = r.getEntry(i, i);
final double[] xi = xData[i];
for (int k = 0; k < xi.length; ++k) {
xi[k] /= rii;
final double lastX = xi[k];
for (int j = i - 1; j >= 0; --j) {
xData[j][k] -= lastX * r.getEntry(j, i);
}
}
}
return new RealMatrixImpl(xData, false);
}
/**
* Check if the decomposed matrix is non-singular.
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular() {
final RealMatrix r = decomposition.getR();
final int p = Math.min(r.getRowDimension(), r.getColumnDimension());
for (int i = 0; i < p; ++i) {
if (r.getEntry(i, i) == 0) {
return false;
}
}
return true;
}
/** Get the pseudo-inverse of the decomposed matrix.
* @return inverse matrix
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealMatrix getInverse()
throws InvalidMatrixException {
final RealMatrix r = decomposition.getR();
final int p = Math.min(r.getRowDimension(), r.getColumnDimension());
return solve(MatrixUtils.createRealIdentityMatrix(p));
}
}

61
src/java/org/apache/commons/math/linear/RealMatrixImpl.java
View File

@ -54,20 +54,15 @@ import org.apache.commons.math.util.MathUtils;
public class RealMatrixImpl implements RealMatrix, Serializable {
/** Serializable version identifier */
private static final long serialVersionUID = 4970229902484487012L;
private static final long serialVersionUID = -391443069570048115L;
/** Entries of the matrix */
protected double data[][];
/** Cached decomposition solver.
/** Cached LU solver.
* @deprecated as of release 2.0, since all methods using this are deprecated
*/
private DecompositionSolver ds;
/** Cached LU decomposition.
* @deprecated as of release 2.0, since all methods using this are deprecated
*/
private LUDecomposition lu;
private LUSolver lu;
/**
* Creates a matrix with no data
@ -89,7 +84,7 @@ public class RealMatrixImpl implements RealMatrix, Serializable {
"row and column dimensions must be postive");
}
data = new double[rowDimension][columnDimension];
ds = null;
lu = null;
}
/**
@ -107,7 +102,7 @@ public class RealMatrixImpl implements RealMatrix, Serializable {
*/
public RealMatrixImpl(double[][] d) {
copyIn(d);
ds = null;
lu = null;
}
/**
@ -147,7 +142,7 @@ public class RealMatrixImpl implements RealMatrix, Serializable {
}
data = d;
}
ds = null;
lu = null;
}
/**
@ -546,7 +541,7 @@ public class RealMatrixImpl implements RealMatrix, Serializable {
System.arraycopy(subMatrix[i], 0, data[row + i], column, nCols);
}
ds = null;
lu = null;
}
@ -631,21 +626,19 @@ public class RealMatrixImpl implements RealMatrix, Serializable {
/** {@inheritDoc} */
@Deprecated
public RealMatrix inverse() throws InvalidMatrixException {
if (ds == null) {
ds = new DecompositionSolver(this);
lu = ds.luDecompose(MathUtils.SAFE_MIN);
if (lu == null) {
lu = new LUSolver(new LUDecompositionImpl(this, MathUtils.SAFE_MIN));
}
return ds.getInverse(lu);
return lu.getInverse();
}
/** {@inheritDoc} */
@Deprecated
public double getDeterminant() throws InvalidMatrixException {
if (ds == null) {
ds = new DecompositionSolver(this);
lu = ds.luDecompose(MathUtils.SAFE_MIN);
if (lu == null) {
lu = new LUSolver(new LUDecompositionImpl(this, MathUtils.SAFE_MIN));
}
return ds.getDeterminant(lu);
return lu.getDeterminant();
}
/** {@inheritDoc} */
@ -656,11 +649,10 @@ public class RealMatrixImpl implements RealMatrix, Serializable {
/** {@inheritDoc} */
@Deprecated
public boolean isSingular() {
if (ds == null) {
ds = new DecompositionSolver(this);
lu = ds.luDecompose(MathUtils.SAFE_MIN);
}
return !ds.isNonSingular(lu);
if (lu == null) {
lu = new LUSolver(new LUDecompositionImpl(this, MathUtils.SAFE_MIN));
}
return !lu.isNonSingular();
}
/** {@inheritDoc} */
@ -792,21 +784,19 @@ public class RealMatrixImpl implements RealMatrix, Serializable {
/** {@inheritDoc} */
@Deprecated
public double[] solve(double[] b) throws IllegalArgumentException, InvalidMatrixException {
if (ds == null) {
ds = new DecompositionSolver(this);
lu = ds.luDecompose(MathUtils.SAFE_MIN);
if (lu == null) {
lu = new LUSolver(new LUDecompositionImpl(this, MathUtils.SAFE_MIN));
}
return ds.solve(b, lu);
return lu.solve(b);
}
/** {@inheritDoc} */
@Deprecated
public RealMatrix solve(RealMatrix b) throws IllegalArgumentException, InvalidMatrixException {
if (ds == null) {
ds = new DecompositionSolver(this);
lu = ds.luDecompose(MathUtils.SAFE_MIN);
if (lu == null) {
lu = new LUSolver(new LUDecompositionImpl(this, MathUtils.SAFE_MIN));
}
return ds.solve(b, lu);
return lu.solve(b);
}
/**
@ -830,9 +820,8 @@ public class RealMatrixImpl implements RealMatrix, Serializable {
*/
@Deprecated
public void luDecompose() throws InvalidMatrixException {
if (ds == null) {
ds = new DecompositionSolver(this);
lu = ds.luDecompose(MathUtils.SAFE_MIN);
if (lu == null) {
lu = new LUSolver(new LUDecompositionImpl(this, MathUtils.SAFE_MIN));
}
}

158
src/java/org/apache/commons/math/linear/SingularValueSolver.java Normal file
View File

@ -0,0 +1,158 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
/**
* Class using singular value decomposition decomposition to solve A &times;
* X = B in least square sense for any matrices A.
* <p>This class solve A &times; X = B in least squares sense: it finds X
* such that ||A &times; X - B|| is minimal.</p>
*
* @version $Revision: 723496 $ $Date: 2008-12-05 00:48:18 +0100 (ven., 05 déc. 2008) $
* @since 2.0
*/
public class SingularValueSolver implements DecompositionSolver {
/** Serializable version identifier. */
private static final long serialVersionUID = -33167987924870528L;
/** Underlying decomposition. */
private final SingularValueDecomposition decomposition;
/**
* Simple constructor.
* @param decomposition decomposition to use
*/
public SingularValueSolver(final SingularValueDecomposition decomposition) {
this.decomposition = decomposition;
}
/** Solve the linear equation A &times; X = B in least square sense.
* <p>The m&times;n matrix A may not be square, the solution X is
* such that ||A &times; X - B|| is minimal.</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 {
final double[] singularValues = decomposition.getSingularValues();
if (b.length != singularValues.length) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] w = decomposition.getUT().operate(b);
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
w[i] /= si;
}
return decomposition.getV().operate(w);
}
/** Solve the linear equation A &times; X = B in least square sense.
* <p>The m&times;n matrix A may not be square, the solution X is
* such that ||A &times; X - B|| is minimal.</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 {
final double[] singularValues = decomposition.getSingularValues();
if (b.getDimension() != singularValues.length) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final RealVector w = decomposition.getUT().operate(b);
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
w.set(i, w.getEntry(i) / si);
}
return decomposition.getV().operate(w);
}
/** Solve the linear equation A &times; X = B in least square sense.
* <p>The m&times;n matrix A may not be square, the solution X is
* such that ||A &times; X - B|| is minimal.</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 {
final double[] singularValues = decomposition.getSingularValues();
if (b.getRowDimension() != singularValues.length) {
throw new IllegalArgumentException("Incorrect row dimension");
}
final RealMatrixImpl w = (RealMatrixImpl) decomposition.getUT().multiply(b);
final double[][] wData = w.getDataRef();
for (int i = 0; i < singularValues.length; ++i) {
final double si = singularValues[i];
if (si == 0) {
throw new SingularMatrixException();
}
final double inv = 1.0 / si;
final double[] wi = wData[i];
for (int j = 0; j < b.getColumnDimension(); ++j) {
wi[j] *= inv;
}
}
return decomposition.getV().multiply(w);
}
/**
* Check if the decomposed matrix is non-singular.
* @return true if the decomposed matrix is non-singular
*/
public boolean isNonSingular() {
return decomposition.getRank() == decomposition.getSingularValues().length;
}
/** Get the pseudo-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();
}
return solve(MatrixUtils.createRealIdentityMatrix(decomposition.getSingularValues().length));
}
}

13
src/java/org/apache/commons/math/stat/regression/GLSMultipleLinearRegression.java
View File

@ -16,7 +16,8 @@
*/
package org.apache.commons.math.stat.regression;
import org.apache.commons.math.linear.DecompositionSolver;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.LUSolver;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealMatrixImpl;
@ -78,8 +79,7 @@ public class GLSMultipleLinearRegression extends AbstractMultipleLinearRegressio
*/
protected RealMatrix getOmegaInverse() {
if (OmegaInverse == null) {
DecompositionSolver solver = new DecompositionSolver(Omega);
OmegaInverse = solver.getInverse(solver.luDecompose());
OmegaInverse = new LUSolver(new LUDecompositionImpl(Omega)).getInverse();
}
return OmegaInverse;
}
@ -95,8 +95,8 @@ public class GLSMultipleLinearRegression extends AbstractMultipleLinearRegressio
RealMatrix OI = getOmegaInverse();
RealMatrix XT = X.transpose();
RealMatrix XTOIX = XT.multiply(OI).multiply(X);
DecompositionSolver solver = new DecompositionSolver(XTOIX);
return solver.getInverse(solver.luDecompose()).multiply(XT).multiply(OI).multiply(Y);
RealMatrix inverse = new LUSolver(new LUDecompositionImpl(XTOIX)).getInverse();
return inverse.multiply(XT).multiply(OI).multiply(Y);
}
/**
@ -109,8 +109,7 @@ public class GLSMultipleLinearRegression extends AbstractMultipleLinearRegressio
protected RealMatrix calculateBetaVariance() {
RealMatrix OI = getOmegaInverse();
RealMatrix XTOIX = X.transpose().multiply(OI).multiply(X);
DecompositionSolver solver = new DecompositionSolver(XTOIX);
return solver.getInverse(solver.luDecompose());
return new LUSolver(new LUDecompositionImpl(XTOIX)).getInverse();
}
/**

5
src/java/org/apache/commons/math/stat/regression/OLSMultipleLinearRegression.java
View File

@ -16,8 +16,8 @@
*/
package org.apache.commons.math.stat.regression;
import org.apache.commons.math.linear.DecompositionSolver;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.LUSolver;
import org.apache.commons.math.linear.QRDecomposition;
import org.apache.commons.math.linear.QRDecompositionImpl;
import org.apache.commons.math.linear.RealMatrix;
@ -109,8 +109,7 @@ public class OLSMultipleLinearRegression extends AbstractMultipleLinearRegressio
*/
protected RealMatrix calculateBetaVariance() {
RealMatrix XTX = X.transpose().multiply(X);
DecompositionSolver solver = new DecompositionSolver(XTX);
return solver.getInverse(solver.luDecompose());
return new LUSolver(new LUDecompositionImpl(XTX)).getInverse();
}

30
src/site/xdoc/userguide/linear.xml
View File

@ -66,8 +66,7 @@ System.out.println(p.getRowDimension()); // 2
System.out.println(p.getColumnDimension()); // 2
// Invert p, using LU decomposition
DecompositionSolver solver = new DecompositionSolver(p);
RealMatrix pInverse = solver.inverse(solver.luDecompose());
RealMatrix pInverse = new LUSolver(new LUDecompositionImpl(p))).getInverse();
</source>
</p>
</subsection>
@ -93,7 +92,7 @@ RealMatrix pInverse = solver.inverse(solver.luDecompose());
</subsection>
<subsection name="3.4 Solving linear systems" href="solve">
<p>
The <code>solve()</code> methods of the <code>DecompositionSolver</code> class
The <code>solve()</code> methods of the <code>DecompositionSolver</code> interface
support solving linear systems of equations of the form AX=B, either in linear sense
or in least square sense. In each case, the <code>RealMatrix</code> represents the
coefficient matrix of the system. For example, to solve the linear system
@ -108,14 +107,13 @@ RealMatrix pInverse = solver.inverse(solver.luDecompose());
RealMatrix coefficients =
new RealMatrixImpl(new double[][] { { 2, 3, -2 }, { -1, 7, 6 }, { 4, -3, -5 } },
false);
DecompositionSolver solver = new DecompositionSolver(coefficients);
LUSolver solver = new LUSolver(new LUDecompositionImpl(coefficients));
</source>
Next create a <code>RealVector</code> array to represent the constant
vector B, select one of the decomposition algorithms (LU, QR ...) and use
<code>solve(vector, decomposition)</code> to solve the system
vector B and use <code>solve(RealVector)</code> to solve the system
<source>
RealVector constants = new RealVectorImpl(new double[] { 1, -2, 1 }, false);
RealVector solution = solver.solve(constants, solver.luDecomposition());
RealVector solution = solver.solve(constants);
</source>
The <code>solution</code> vector will contain values for x
(<code>solution.getEntry(0)</code>), y (<code>solution.getEntry(1)</code>),
@ -126,11 +124,10 @@ RealVector solution = solver.solve(constants, solver.luDecomposition());
for X is such that the residual AX-B has minimal norm. If an exact solution
exist (i.e. if for some X the residual AX-B is exactly 0), then this exact
solution is also the solution in least square sense. Some solvers like
<code>LUDecomposition</code> can only find the solution for
square matrices and when the solution is an exact linear solution. Other solvers
like <code>QRDecomposition</code> are more versatile and can also find solutions
with non-square matrix A or when no exact solution exist (i.e. when the minimal
value for AX-B norm is non-null).
<code>LUSolver</code> can only find the solution for square matrices and when
the solution is an exact linear solution. Other solvers like <code>QRDecomposition</code>
are more versatile and can also find solutions with non-square matrix A or when
no exact solution exist (i.e. when the minimal value for AX-B norm is non-null).
</p>
<p>
If the coefficient matrix is singular or doesn't fulfill the requirements of
@ -145,12 +142,15 @@ RealVector solution = solver.solve(constants, solver.luDecomposition());
<p>
It is possible to solve multiple systems with the same coefficient matrix
in one method call. To do this, create a matrix whose column vectors correspond
to the constant vectors for the systems to be solved and use
<code>solve(matrix, decomposition),</code> which returns a matrix with column
vectors representing the solutions.
to the constant vectors for the systems to be solved and use <code>solve(RealMatrix),</code>
which returns a matrix with column vectors representing the solutions.
</p>
</subsection>
<subsection name="3.5 Eigenvalues/eigenvectors and singular values/singular vectors" href="eigen">
<p>
Decomposition algorithms may be used for themselves and not only for linear system solving.
This is of prime interest with eigen decomposition and singular value decomposition.
</p>
<p>
The <code>getEigenvalue()</code>, <code>getEigenvalues()</code>, <code>getEigenVector()</code>,
<code>getV()</code>, <code>getD()</code> and <code>getVT()</code> methods of the

150
src/test/org/apache/commons/math/linear/EigenDecompositionImplTest.java
View File

@ -46,7 +46,7 @@ public class EigenDecompositionImplTest extends TestCase {
new RealMatrixImpl(new double[][] {
{ 1.5 }
}, false);
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(1.5, ed.getEigenvalue(0), 1.0e-15);
}
@ -56,7 +56,7 @@ public class EigenDecompositionImplTest extends TestCase {
{ 59.0, 12.0 },
{ Double.NaN, 66.0 }
}, false);
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(75.0, ed.getEigenvalue(0), 1.0e-15);
assertEquals(50.0, ed.getEigenvalue(1), 1.0e-15);
}
@ -68,7 +68,7 @@ public class EigenDecompositionImplTest extends TestCase {
{ Double.NaN, 8693.0, 7920.0 },
{ Double.NaN, Double.NaN, 17300.0 }
}, false);
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(50000.0, ed.getEigenvalue(0), 3.0e-11);
assertEquals(12500.0, ed.getEigenvalue(1), 3.0e-11);
assertEquals( 3125.0, ed.getEigenvalue(2), 3.0e-11);
@ -82,7 +82,7 @@ public class EigenDecompositionImplTest extends TestCase {
{ Double.NaN, Double.NaN, 0.164, -0.048 },
{ Double.NaN, Double.NaN, Double.NaN, 0.136 }
}, false);
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(1.0, ed.getEigenvalue(0), 1.0e-15);
assertEquals(0.4, ed.getEigenvalue(1), 1.0e-15);
assertEquals(0.2, ed.getEigenvalue(2), 1.0e-15);
@ -97,7 +97,7 @@ public class EigenDecompositionImplTest extends TestCase {
{ 0.1152, -0.2304, 0.3088, -0.1344 },
{ -0.2976, 0.1152, -0.1344, 0.3872 }
}, false);
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(1.0, ed.getEigenvalue(0), 1.0e-15);
assertEquals(0.4, ed.getEigenvalue(1), 1.0e-15);
assertEquals(0.2, ed.getEigenvalue(2), 1.0e-15);
@ -133,7 +133,7 @@ public class EigenDecompositionImplTest extends TestCase {
/** test dimensions */
public void testDimensions() {
final int m = matrix.getRowDimension();
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
assertEquals(m, ed.getV().getRowDimension());
assertEquals(m, ed.getV().getColumnDimension());
assertEquals(m, ed.getD().getColumnDimension());
@ -144,7 +144,7 @@ public class EigenDecompositionImplTest extends TestCase {
/** test eigenvalues */
public void testEigenvalues() {
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
double[] eigenValues = ed.getEigenvalues();
assertEquals(refValues.length, eigenValues.length);
for (int i = 0; i < refValues.length; ++i) {
@ -161,7 +161,7 @@ public class EigenDecompositionImplTest extends TestCase {
}
Arrays.sort(bigValues);
EigenDecomposition ed =
new DecompositionSolver(createTestMatrix(r, bigValues)).eigenDecompose();
new EigenDecompositionImpl(createTestMatrix(r, bigValues), MathUtils.SAFE_MIN);
double[] eigenValues = ed.getEigenvalues();
assertEquals(bigValues.length, eigenValues.length);
for (int i = 0; i < bigValues.length; ++i) {
@ -171,7 +171,7 @@ public class EigenDecompositionImplTest extends TestCase {
/** test eigenvectors */
public void testEigenvectors() {
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
for (int i = 0; i < matrix.getRowDimension(); ++i) {
double lambda = ed.getEigenvalue(i);
RealVector v = ed.getEigenvector(i);
@ -182,7 +182,7 @@ public class EigenDecompositionImplTest extends TestCase {
/** test A = VDVt */
public void testAEqualVDVt() {
EigenDecomposition ed = new DecompositionSolver(matrix).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN);
RealMatrix v = ed.getV();
RealMatrix d = ed.getD();
RealMatrix vT = ed.getVT();
@ -192,143 +192,23 @@ public class EigenDecompositionImplTest extends TestCase {
/** test that V is orthogonal */
public void testVOrthogonal() {
RealMatrix v = new DecompositionSolver(matrix).eigenDecompose().getV();
RealMatrix v = new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN).getV();
RealMatrix vTv = v.transpose().multiply(v);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(vTv.getRowDimension());
assertEquals(0, vTv.subtract(id).getNorm(), 2.0e-13);
}
/** test non invertible matrix */
public void testNonInvertible() {
Random r = new Random(9994100315209l);
DecompositionSolver ds =
new DecompositionSolver(createTestMatrix(r, new double[] { 1.0, 0.0, -1.0, -2.0, -3.0 }));
EigenDecomposition ed = ds.eigenDecompose();
assertFalse(ds.isNonSingular(ed));
try {
ds.getInverse(ed);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test invertible matrix */
public void testInvertible() {
Random r = new Random(9994100315209l);
RealMatrix m =
createTestMatrix(r, new double[] { 1.0, 0.5, -1.0, -2.0, -3.0 });
DecompositionSolver ds = new DecompositionSolver(m);
EigenDecomposition ed = ds.eigenDecompose();
assertTrue(ds.isNonSingular(ed));
RealMatrix inverse = ds.getInverse(ed);
RealMatrix error =
m.multiply(inverse).subtract(MatrixUtils.createRealIdentityMatrix(m.getRowDimension()));
assertEquals(0, error.getNorm(), 4.0e-15);
}
/** test diagonal matrix */
public void testDiagonal() {
double[] diagonal = new double[] { -3.0, -2.0, 2.0, 5.0 };
RealMatrix m = createDiagonalMatrix(diagonal, diagonal.length, diagonal.length);
EigenDecomposition ed = new DecompositionSolver(m).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(m, MathUtils.SAFE_MIN);
assertEquals(diagonal[0], ed.getEigenvalue(3), 2.0e-15);
assertEquals(diagonal[1], ed.getEigenvalue(2), 2.0e-15);
assertEquals(diagonal[2], ed.getEigenvalue(1), 2.0e-15);
assertEquals(diagonal[3], ed.getEigenvalue(0), 2.0e-15);
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
DecompositionSolver ds = new DecompositionSolver(matrix);
EigenDecomposition ed = ds.eigenDecompose();
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
ds.solve(b, ed);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ds.solve(b.getColumn(0), ed);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ds.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)), ed);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
RealMatrix m = new RealMatrixImpl(new double[][] {
{ 91, 5, 29, 32, 40, 14 },
{ 5, 34, -1, 0, 2, -1 },
{ 29, -1, 12, 9, 21, 8 },
{ 32, 0, 9, 14, 9, 0 },
{ 40, 2, 21, 9, 51, 19 },
{ 14, -1, 8, 0, 19, 14 }
});
DecompositionSolver ds = new DecompositionSolver(m);
EigenDecomposition ed = ds.eigenDecompose();
assertEquals(184041, ds.getDeterminant(ed), 2.0e-8);
RealMatrix b = new RealMatrixImpl(new double[][] {
{ 1561, 269, 188 },
{ 69, -21, 70 },
{ 739, 108, 63 },
{ 324, 86, 59 },
{ 1624, 194, 107 },
{ 796, 69, 36 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ 1, 2, 1 },
{ 2, -1, 2 },
{ 4, 2, 3 },
{ 8, -1, 0 },
{ 16, 2, 0 },
{ 32, -1, 0 }
});
// using RealMatrix
assertEquals(0, ds.solve(b, ed).subtract(xRef).getNorm(), 2.0e-12);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(ds.solve(b.getColumn(i), ed)).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
// using RealMatrixImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
ds.solve(b.getColumnVector(i), ed).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
// using RealMatrix with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
ds.solve(v, ed).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
}
/**
* Matrix with eigenvalues {8, -1, -1}
*/
@ -338,7 +218,7 @@ public class EigenDecompositionImplTest extends TestCase {
{2, 0, 2},
{4, 2, 3}
});
EigenDecomposition ed = new DecompositionSolver(repeated).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(repeated, MathUtils.SAFE_MIN);
checkEigenValues((new double[] {8, -1, -1}), ed, 1E-12);
checkEigenVector((new double[] {2, 1, 2}), ed, 1E-12);
}
@ -352,7 +232,7 @@ public class EigenDecompositionImplTest extends TestCase {
{1, 3, -4},
{-4, -4, 8}
});
EigenDecomposition ed = new DecompositionSolver(distinct).eigenDecompose();
EigenDecomposition ed = new EigenDecompositionImpl(distinct, MathUtils.SAFE_MIN);
checkEigenValues((new double[] {2, 0, 12}), ed, 1E-12);
checkEigenVector((new double[] {1, -1, 0}), ed, 1E-12);
checkEigenVector((new double[] {1, 1, 1}), ed, 1E-12);
@ -442,7 +322,7 @@ public class EigenDecompositionImplTest extends TestCase {
matrix = null;
}
private RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
static RealMatrix createTestMatrix(final Random r, final double[] eigenValues) {
final int n = eigenValues.length;
final RealMatrix v = createOrthogonalMatrix(r, n);
final RealMatrix d = createDiagonalMatrix(eigenValues, n, n);

172
src/test/org/apache/commons/math/linear/EigenSolverTest.java Normal file
View File

@ -0,0 +1,172 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import java.util.Random;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
import org.apache.commons.math.util.MathUtils;
public class EigenSolverTest extends TestCase {
private double[] refValues;
private RealMatrix matrix;
public EigenSolverTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(EigenSolverTest.class);
suite.setName("EigenSolver Tests");
return suite;
}
/** test non invertible matrix */
public void testNonInvertible() {
Random r = new Random(9994100315209l);
RealMatrix m =
EigenDecompositionImplTest.createTestMatrix(r, new double[] { 1.0, 0.0, -1.0, -2.0, -3.0 });
EigenSolver es = new EigenSolver(new EigenDecompositionImpl(m, MathUtils.SAFE_MIN));
assertFalse(es.isNonSingular());
try {
es.getInverse();
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test invertible matrix */
public void testInvertible() {
Random r = new Random(9994100315209l);
RealMatrix m =
EigenDecompositionImplTest.createTestMatrix(r, new double[] { 1.0, 0.5, -1.0, -2.0, -3.0 });
EigenSolver es = new EigenSolver(new EigenDecompositionImpl(m, MathUtils.SAFE_MIN));
assertTrue(es.isNonSingular());
RealMatrix inverse = es.getInverse();
RealMatrix error =
m.multiply(inverse).subtract(MatrixUtils.createRealIdentityMatrix(m.getRowDimension()));
assertEquals(0, error.getNorm(), 4.0e-15);
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
EigenSolver es = new EigenSolver(new EigenDecompositionImpl(matrix, MathUtils.SAFE_MIN));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
es.solve(b);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
es.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
es.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
RealMatrix m = new RealMatrixImpl(new double[][] {
{ 91, 5, 29, 32, 40, 14 },
{ 5, 34, -1, 0, 2, -1 },
{ 29, -1, 12, 9, 21, 8 },
{ 32, 0, 9, 14, 9, 0 },
{ 40, 2, 21, 9, 51, 19 },
{ 14, -1, 8, 0, 19, 14 }
});
EigenSolver es = new EigenSolver(new EigenDecompositionImpl(m, MathUtils.SAFE_MIN));
assertEquals(184041, es.getDeterminant(), 2.0e-8);
RealMatrix b = new RealMatrixImpl(new double[][] {
{ 1561, 269, 188 },
{ 69, -21, 70 },
{ 739, 108, 63 },
{ 324, 86, 59 },
{ 1624, 194, 107 },
{ 796, 69, 36 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ 1, 2, 1 },
{ 2, -1, 2 },
{ 4, 2, 3 },
{ 8, -1, 0 },
{ 16, 2, 0 },
{ 32, -1, 0 }
});
// using RealMatrix
assertEquals(0, es.solve(b).subtract(xRef).getNorm(), 2.0e-12);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(es.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
// using RealMatrixImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
es.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
// using RealMatrix with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
es.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
2.0e-11);
}
}
public void setUp() {
refValues = new double[] {
2.003, 2.002, 2.001, 1.001, 1.000, 0.001
};
matrix = EigenDecompositionImplTest.createTestMatrix(new Random(35992629946426l), refValues);
}
public void tearDown() {
refValues = null;
matrix = null;
}
}

137
src/test/org/apache/commons/math/linear/LUDecompositionImplTest.java
View File

@ -88,18 +88,6 @@ public class LUDecompositionImplTest extends TestCase {
}
}
/** test threshold impact */
public void testThreshold() {
final RealMatrix matrix = new RealMatrixImpl(new double[][] {
{ 1.0, 2.0, 3.0},
{ 2.0, 5.0, 3.0},
{ 4.000001, 9.0, 9.0}
}, false);
DecompositionSolver solver = new DecompositionSolver(matrix);
assertFalse(solver.isNonSingular(solver.luDecompose(1.0e-5)));
assertTrue(solver.isNonSingular(solver.luDecompose(1.0e-10)));
}
/** test PA = LU */
public void testPAEqualLU() {
RealMatrix matrix = new RealMatrixImpl(testData, false);
@ -226,118 +214,6 @@ public class LUDecompositionImplTest extends TestCase {
assertTrue(lu.isSingular());
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
DecompositionSolver solver =
new DecompositionSolver(new RealMatrixImpl(testData, false));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
solver.solve(b, solver.luDecompose());
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0), solver.luDecompose());
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)),
solver.luDecompose());
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve singularity errors */
public void testSolveSingularityErrors() {
DecompositionSolver solver =
new DecompositionSolver(new RealMatrixImpl(singular, false));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
solver.solve(b, solver.luDecompose());
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0), solver.luDecompose());
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0), solver.luDecompose());
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)),
solver.luDecompose());
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
DecompositionSolver solver =
new DecompositionSolver(new RealMatrixImpl(testData, false));
LUDecomposition lu = solver.luDecompose();
RealMatrix b = new RealMatrixImpl(new double[][] {
{ 1, 0 }, { 2, -5 }, { 3, 1 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ 19, -71 }, { -6, 22 }, { -2, 9 }
});
// using RealMatrix
assertEquals(0, solver.solve(b, lu).subtract(xRef).getNorm(), 1.0e-13);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(solver.solve(b.getColumn(i), lu)).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVectorImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
solver.solve(b.getColumnVector(i), lu).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVector with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
solver.solve(v, lu).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
}
/** test matrices values */
public void testMatricesValues1() {
LUDecomposition lu =
@ -418,17 +294,4 @@ public class LUDecompositionImplTest extends TestCase {
}
/** test determinant */
public void testDeterminant() {
assertEquals( -1, getDeterminant(new RealMatrixImpl(testData, false)), 1.0e-15);
assertEquals(-10, getDeterminant(new RealMatrixImpl(luData, false)), 1.0e-14);
assertEquals( 0, getDeterminant(new RealMatrixImpl(singular, false)), 1.0e-17);
assertEquals( 0, getDeterminant(new RealMatrixImpl(bigSingular, false)), 1.0e-10);
}
private double getDeterminant(RealMatrix m) {
DecompositionSolver ds = new DecompositionSolver(m);
return ds.getDeterminant(ds.luDecompose());
}
}

201
src/test/org/apache/commons/math/linear/LUSolverTest.java Normal file
View File

@ -0,0 +1,201 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
public class LUSolverTest extends TestCase {
private double[][] testData = {
{ 1.0, 2.0, 3.0},
{ 2.0, 5.0, 3.0},
{ 1.0, 0.0, 8.0}
};
private double[][] luData = {
{ 2.0, 3.0, 3.0 },
{ 0.0, 5.0, 7.0 },
{ 6.0, 9.0, 8.0 }
};
// singular matrices
private double[][] singular = {
{ 2.0, 3.0 },
{ 2.0, 3.0 }
};
private double[][] bigSingular = {
{ 1.0, 2.0, 3.0, 4.0 },
{ 2.0, 5.0, 3.0, 4.0 },
{ 7.0, 3.0, 256.0, 1930.0 },
{ 3.0, 7.0, 6.0, 8.0 }
}; // 4th row = 1st + 2nd
public LUSolverTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(LUSolverTest.class);
suite.setName("LUSolver Tests");
return suite;
}
/** test threshold impact */
public void testThreshold() {
final RealMatrix matrix = new RealMatrixImpl(new double[][] {
{ 1.0, 2.0, 3.0},
{ 2.0, 5.0, 3.0},
{ 4.000001, 9.0, 9.0}
}, false);
assertFalse(new LUSolver(new LUDecompositionImpl(matrix, 1.0e-5)).isNonSingular());
assertTrue(new LUSolver(new LUDecompositionImpl(matrix, 1.0e-10)).isNonSingular());
}
/** test singular */
public void testSingular() {
LUSolver lu =
new LUSolver(new LUDecompositionImpl(new RealMatrixImpl(testData, false)));
assertTrue(lu.isNonSingular());
lu = new LUSolver(new LUDecompositionImpl(new RealMatrixImpl(singular, false)));
assertFalse(lu.isNonSingular());
lu = new LUSolver(new LUDecompositionImpl(new RealMatrixImpl(bigSingular, false)));
assertFalse(lu.isNonSingular());
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
LUSolver solver =
new LUSolver(new LUDecompositionImpl(new RealMatrixImpl(testData, false)));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve singularity errors */
public void testSolveSingularityErrors() {
LUSolver solver =
new LUSolver(new LUDecompositionImpl(new RealMatrixImpl(singular, false)));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
LUSolver solver =
new LUSolver(new LUDecompositionImpl(new RealMatrixImpl(testData, false)));
RealMatrix b = new RealMatrixImpl(new double[][] {
{ 1, 0 }, { 2, -5 }, { 3, 1 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ 19, -71 }, { -6, 22 }, { -2, 9 }
});
// using RealMatrix
assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVectorImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVector with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
}
/** test determinant */
public void testDeterminant() {
assertEquals( -1, getDeterminant(new RealMatrixImpl(testData, false)), 1.0e-15);
assertEquals(-10, getDeterminant(new RealMatrixImpl(luData, false)), 1.0e-14);
assertEquals( 0, getDeterminant(new RealMatrixImpl(singular, false)), 1.0e-17);
assertEquals( 0, getDeterminant(new RealMatrixImpl(bigSingular, false)), 1.0e-10);
}
private double getDeterminant(RealMatrix m) {
return new LUSolver(new LUDecompositionImpl(m)).getDeterminant();
}
}

124
src/test/org/apache/commons/math/linear/QRDecompositionImplTest.java
View File

@ -196,130 +196,10 @@ public class QRDecompositionImplTest extends TestCase {
}
/** test rank */
public void testRank() {
DecompositionSolver ds =
new DecompositionSolver(new RealMatrixImpl(testData3x3NonSingular, false));
assertTrue(ds.isNonSingular(ds.qrDecompose()));
ds = new DecompositionSolver(new RealMatrixImpl(testData3x3Singular, false));
assertFalse(ds.isNonSingular(ds.qrDecompose()));
ds = new DecompositionSolver(new RealMatrixImpl(testData3x4, false));
assertTrue(ds.isNonSingular(ds.qrDecompose()));
ds = new DecompositionSolver(new RealMatrixImpl(testData4x3, false));
assertTrue(ds.isNonSingular(ds.qrDecompose()));
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
DecompositionSolver solver =
new DecompositionSolver(new RealMatrixImpl(testData3x3NonSingular, false));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
solver.solve(b, solver.qrDecompose());
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0), solver.qrDecompose());
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0), solver.qrDecompose());
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve rank errors */
public void testSolveRankErrors() {
DecompositionSolver solver =
new DecompositionSolver(new RealMatrixImpl(testData3x3Singular, false));
RealMatrix b = new RealMatrixImpl(new double[3][2]);
try {
solver.solve(b, solver.qrDecompose());
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0), solver.qrDecompose());
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0), solver.qrDecompose());
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
DecompositionSolver ds =
new DecompositionSolver(new RealMatrixImpl(testData3x3NonSingular, false));
QRDecomposition qr = ds.qrDecompose();
RealMatrix b = new RealMatrixImpl(new double[][] {
{ -102, 12250 }, { 544, 24500 }, { 167, -36750 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ 1, 2515 }, { 2, 422 }, { -3, 898 }
});
// using RealMatrix
assertEquals(0, ds.solve(b, qr).subtract(xRef).getNorm(), 1.0e-13);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(ds.solve(b.getColumn(i), qr)).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVectorImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
ds.solve(b.getColumnVector(i), qr).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVector with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
ds.solve(v, qr).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
}
/** test matrices values */
public void testMatricesValues() {
DecompositionSolver ds =
new DecompositionSolver(new RealMatrixImpl(testData3x3NonSingular, false));
QRDecomposition qr = ds.qrDecompose();
QRDecomposition qr =
new QRDecompositionImpl(new RealMatrixImpl(testData3x3NonSingular, false));
RealMatrix qRef = new RealMatrixImpl(new double[][] {
{ -12.0 / 14.0, 69.0 / 175.0, -58.0 / 175.0 },
{ -6.0 / 14.0, -158.0 / 175.0, 6.0 / 175.0 },

174
src/test/org/apache/commons/math/linear/QRSolverTest.java Normal file
View File

@ -0,0 +1,174 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
public class QRSolverTest extends TestCase {
double[][] testData3x3NonSingular = {
{ 12, -51, 4 },
{ 6, 167, -68 },
{ -4, 24, -41 }, };
double[][] testData3x3Singular = {
{ 1, 4, 7, },
{ 2, 5, 8, },
{ 3, 6, 9, }, };
double[][] testData3x4 = {
{ 12, -51, 4, 1 },
{ 6, 167, -68, 2 },
{ -4, 24, -41, 3 }, };
double[][] testData4x3 = {
{ 12, -51, 4, },
{ 6, 167, -68, },
{ -4, 24, -41, },
{ -5, 34, 7, }, };
public QRSolverTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(QRSolverTest.class);
suite.setName("QRSolver Tests");
return suite;
}
/** test rank */
public void testRank() {
QRSolver solver =
new QRSolver(new QRDecompositionImpl(new RealMatrixImpl(testData3x3NonSingular, false)));
assertTrue(solver.isNonSingular());
solver = new QRSolver(new QRDecompositionImpl(new RealMatrixImpl(testData3x3Singular, false)));
assertFalse(solver.isNonSingular());
solver = new QRSolver(new QRDecompositionImpl(new RealMatrixImpl(testData3x4, false)));
assertTrue(solver.isNonSingular());
solver = new QRSolver(new QRDecompositionImpl(new RealMatrixImpl(testData4x3, false)));
assertTrue(solver.isNonSingular());
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
QRSolver solver =
new QRSolver(new QRDecompositionImpl(new RealMatrixImpl(testData3x3NonSingular, false)));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve rank errors */
public void testSolveRankErrors() {
QRSolver solver =
new QRSolver(new QRDecompositionImpl(new RealMatrixImpl(testData3x3Singular, false)));
RealMatrix b = new RealMatrixImpl(new double[3][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
QRSolver solver =
new QRSolver(new QRDecompositionImpl(new RealMatrixImpl(testData3x3NonSingular, false)));
RealMatrix b = new RealMatrixImpl(new double[][] {
{ -102, 12250 }, { 544, 24500 }, { 167, -36750 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ 1, 2515 }, { 2, 422 }, { -3, 898 }
});
// using RealMatrix
assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVectorImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealVector with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
}
}

12
src/test/org/apache/commons/math/linear/RealMatrixImplTest.java
View File

@ -261,10 +261,8 @@ public final class RealMatrixImplTest extends TestCase {
/** test transpose */
public void testTranspose() {
RealMatrix m = new RealMatrixImpl(testData);
DecompositionSolver ds1 = new DecompositionSolver(m);
RealMatrix mIT = ds1.getInverse(ds1.luDecompose()).transpose();
DecompositionSolver ds2 = new DecompositionSolver(m.transpose());
RealMatrix mTI = ds2.getInverse(ds2.luDecompose());
RealMatrix mIT = new LUSolver(new LUDecompositionImpl(m)).getInverse().transpose();
RealMatrix mTI = new LUSolver(new LUDecompositionImpl(m.transpose())).getInverse();
assertClose("inverse-transpose", mIT, mTI, normTolerance);
m = new RealMatrixImpl(testData2);
RealMatrix mt = new RealMatrixImpl(testData2T);
@ -354,8 +352,7 @@ public final class RealMatrixImplTest extends TestCase {
assertEquals(2, p.getRowDimension());
assertEquals(2, p.getColumnDimension());
// Invert p
DecompositionSolver ds1 = new DecompositionSolver(p);
RealMatrix pInverse = ds1.getInverse(ds1.luDecompose());
RealMatrix pInverse = new LUSolver(new LUDecompositionImpl(p)).getInverse();
assertEquals(2, pInverse.getRowDimension());
assertEquals(2, pInverse.getColumnDimension());
@ -363,8 +360,7 @@ public final class RealMatrixImplTest extends TestCase {
double[][] coefficientsData = {{2, 3, -2}, {-1, 7, 6}, {4, -3, -5}};
RealMatrix coefficients = new RealMatrixImpl(coefficientsData);
double[] constants = {1, -2, 1};
DecompositionSolver ds2 = new DecompositionSolver(coefficients);
double[] solution = ds2.solve(constants, ds2.luDecompose());
double[] solution = new LUSolver(new LUDecompositionImpl(coefficients)).solve(constants);
assertEquals(2 * solution[0] + 3 * solution[1] -2 * solution[2], constants[0], 1E-12);
assertEquals(-1 * solution[0] + 7 * solution[1] + 6 * solution[2], constants[1], 1E-12);
assertEquals(4 * solution[0] - 3 * solution[1] -5 * solution[2], constants[2], 1E-12);

118
src/test/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java
View File

@ -129,122 +129,6 @@ public class SingularValueDecompositionImplTest extends TestCase {
assertEquals(0, mTm.subtract(id).getNorm(), normTolerance);
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
DecompositionSolver ds =
new DecompositionSolver(new RealMatrixImpl(testSquare, false));
SingularValueDecomposition svd = ds.singularDecompose();
RealMatrix b = new RealMatrixImpl(new double[3][2]);
try {
ds.solve(b, svd);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ds.solve(b.getColumn(0), svd);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ds.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)), svd);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve singularity errors */
public void testSolveSingularityErrors() {
DecompositionSolver ds =
new DecompositionSolver(new RealMatrixImpl(new double[][] {
{ 1.0, 0.0 },
{ 0.0, 0.0 }
}, false));
SingularValueDecomposition svd = ds.singularDecompose();
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
ds.solve(b, svd);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ds.solve(b.getColumn(0), svd);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ds.solve(b.getColumnVector(0), svd);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
ds.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)), svd);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
DecompositionSolver ds =
new DecompositionSolver(new RealMatrixImpl(testSquare, false));
SingularValueDecomposition svd = ds.singularDecompose();
RealMatrix b = new RealMatrixImpl(new double[][] {
{ 1, 2, 3 }, { 0, -5, 1 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ -8.0 / 25.0, -263.0 / 75.0, -29.0 / 75.0 },
{ 19.0 / 25.0, 78.0 / 25.0, 49.0 / 25.0 }
});
// using RealMatrix
assertEquals(0, ds.solve(b, svd).subtract(xRef).getNorm(), normTolerance);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(ds.solve(b.getColumn(i), svd)).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealMatrixImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
ds.solve(b.getColumnVector(i), svd).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealMatrix with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
ds.solve(v, svd).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
}
/** test matrices values */
public void testMatricesValues1() {
SingularValueDecomposition svd =
@ -299,7 +183,7 @@ public class SingularValueDecompositionImplTest extends TestCase {
// check values against known references
SingularValueDecomposition svd =
new DecompositionSolver(new RealMatrixImpl(testNonSquare, false)).singularDecompose();
new SingularValueDecompositionImpl(new RealMatrixImpl(testNonSquare, false));
RealMatrix u = svd.getU();
assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
RealMatrix s = svd.getS();

164
src/test/org/apache/commons/math/linear/SingularValueSolverTest.java Normal file
View File

@ -0,0 +1,164 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.linear;
import junit.framework.Test;
import junit.framework.TestCase;
import junit.framework.TestSuite;
public class SingularValueSolverTest extends TestCase {
private double[][] testSquare = {
{ 24.0 / 25.0, 43.0 / 25.0 },
{ 57.0 / 25.0, 24.0 / 25.0 }
};
private static final double normTolerance = 10e-14;
public SingularValueSolverTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(SingularValueSolverTest.class);
suite.setName("SingularValueSolver Tests");
return suite;
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
SingularValueSolver solver =
new SingularValueSolver(new SingularValueDecompositionImpl(new RealMatrixImpl(testSquare, false)));
RealMatrix b = new RealMatrixImpl(new double[3][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve singularity errors */
public void testSolveSingularityErrors() {
RealMatrix m =
new RealMatrixImpl(new double[][] {
{ 1.0, 0.0 },
{ 0.0, 0.0 }
}, false);
SingularValueSolver solver = new SingularValueSolver(new SingularValueDecompositionImpl(m));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
solver.solve(new RealVectorImplTest.RealVectorTestImpl(b.getColumn(0)));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
}
/** test solve */
public void testSolve() {
SingularValueSolver solver =
new SingularValueSolver(new SingularValueDecompositionImpl(new RealMatrixImpl(testSquare, false)));
RealMatrix b = new RealMatrixImpl(new double[][] {
{ 1, 2, 3 }, { 0, -5, 1 }
});
RealMatrix xRef = new RealMatrixImpl(new double[][] {
{ -8.0 / 25.0, -263.0 / 75.0, -29.0 / 75.0 },
{ 19.0 / 25.0, 78.0 / 25.0, 49.0 / 25.0 }
});
// using RealMatrix
assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), normTolerance);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealMatrixImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealMatrix with an alternate implementation
for (int i = 0; i < b.getColumnDimension(); ++i) {
RealVectorImplTest.RealVectorTestImpl v =
new RealVectorImplTest.RealVectorTestImpl(b.getColumn(i));
assertEquals(0,
solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
}
/** test condition number */
public void testConditionNumber() {
SingularValueDecompositionImpl svd =
new SingularValueDecompositionImpl(new RealMatrixImpl(testSquare, false));
assertEquals(3.0, svd.getConditionNumber(), 1.0e-15);
}
}