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

added (at least!) a first implementation of singular value decomposition. 

JIRA: MATH-157, MATH-220, MATH-233

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@723361 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2008-12-04 15:46:44 +00:00
parent f1fa39e47a
commit 2655bad9b9
4 changed files with 810 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

77
src/java/org/apache/commons/math/linear/SingularValueDecomposition.java
View File

@ -17,7 +17,6 @@
package org.apache.commons.math.linear;
import org.apache.commons.math.ConvergenceException;
/**
* An interface to classes that implement an algorithm to calculate the
@ -34,8 +33,6 @@ import org.apache.commons.math.ConvergenceException;
* <li><code>solve</code> methods have been added (in the superinterface),</li>
* <li>a {@link DecompositionSolver#decompose(RealMatrix) decompose(RealMatrix)}
* method has been added (in the superinterface),</li>
* <li>a {@link #decompose(RealMatrix, int) decompose(RealMatrix), int)} method
* has been added,</li>
* <li>a {@link DecompositionSolver#isNonSingular() isNonSingular} method has
* been added (in the superinterface),</li>
* <li>a {@link DecompositionSolver#getInverse() getInverse} method has been
@ -54,42 +51,43 @@ import org.apache.commons.math.ConvergenceException;
*/
public interface SingularValueDecomposition extends DecompositionSolver {
/**
* Decompose a matrix to find its largest singular values.
* @param matrix matrix to decompose
* @param maxSingularValues maximal number of singular values to compute
* @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
* if algorithm fails to converge
*/
void decompose(RealMatrix matrix, int maxSingularValues)
throws InvalidMatrixException;
/**
* Returns the matrix U of the decomposition.
* <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
* @return the U matrix
* @exception IllegalStateException if neither {@link
* DecompositionSolver#decompose(RealMatrix) decompose} nor {@link
* #decompose(RealMatrix, int)} have not been called
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
* @see #getUT()
*/
RealMatrix getU() throws IllegalStateException;
/**
* Returns the transpose of the matrix U of the decomposition.
* <p>U is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
* @return the U matrix (or null if decomposed matrix is singular)
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
* @see #getU()
*/
RealMatrix getUT() throws IllegalStateException;
/**
* Returns the diagonal matrix &Sigma; of the decomposition.
* <p>&Sigma; is a diagonal matrix.</p>
* <p>&Sigma; is a diagonal matrix. The singular values are provided in
* non-increasing order, for compatibility with Jama.</p>
* @return the &Sigma; matrix
* @exception IllegalStateException if neither {@link
* DecompositionSolver#decompose(RealMatrix) decompose} nor {@link
* #decompose(RealMatrix, int)} have not been called
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
*/
RealMatrix getS() throws IllegalStateException;
/**
* Returns the diagonal elements of the matrix &Sigma; of the decomposition.
* Returns the diagonal elements of the matrix &Sigma; of the decomposition.
* <p>The singular values are provided in non-increasing order, for
* compatibility with Jama.</p>
* @return the diagonal elements of the &Sigma; matrix
* @exception IllegalStateException if neither {@link
* DecompositionSolver#decompose(RealMatrix) decompose} nor {@link
* #decompose(RealMatrix, int)} have not been called
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
*/
double[] getSingularValues() throws IllegalStateException;
@ -97,30 +95,38 @@ public interface SingularValueDecomposition extends DecompositionSolver {
* Returns the matrix V of the decomposition.
* <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
* @return the V matrix (or null if decomposed matrix is singular)
* @exception IllegalStateException if neither {@link
* DecompositionSolver#decompose(RealMatrix) decompose} nor {@link
* #decompose(RealMatrix, int)} have not been called
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
* @see #getVT()
*/
RealMatrix getV() throws IllegalStateException;
/**
* Returns the transpose of the matrix V of the decomposition.
* <p>V is an orthogonal matrix, i.e. its transpose is also its inverse.</p>
* @return the V matrix (or null if decomposed matrix is singular)
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
* @see #getV()
*/
RealMatrix getVT() throws IllegalStateException;
/**
* Returns the L<sub>2</sub> norm of the matrix.
* <p>The L<sub>2</sub> norm is max(|A &times; u|<sub>2</sub> /
* |u|<sub>2</sub>), where |.|<sub>2</sub> denotes the vectorial 2-norm
* (i.e. the traditional euclidian norm).</p>
* @return norm
* @exception IllegalStateException if neither {@link
* DecompositionSolver#decompose(RealMatrix) decompose} nor {@link
* #decompose(RealMatrix, int)} have not been called
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
*/
double getNorm() throws IllegalStateException;
/**
* Return the condition number of the matrix.
* @return condition number of the matrix
* @exception IllegalStateException if neither {@link
* DecompositionSolver#decompose(RealMatrix) decompose} nor {@link
* #decompose(RealMatrix, int)} have not been called
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
*/
double getConditionNumber() throws IllegalStateException;
@ -131,9 +137,8 @@ public interface SingularValueDecomposition extends DecompositionSolver {
* terms is max(m,n) &times; ulp(s<sub>1</sub>) where ulp(s<sub>1</sub>)
* is the least significant bit of the largest singular value.</p>
* @return effective numerical matrix rank
* @exception IllegalStateException if neither {@link
* DecompositionSolver#decompose(RealMatrix) decompose} nor {@link
* #decompose(RealMatrix, int)} have not been called
* @exception IllegalStateException if {@link
* DecompositionSolver#decompose(RealMatrix) decompose} has not been called
*/
int getRank() throws IllegalStateException;

435
src/java/org/apache/commons/math/linear/SingularValueDecompositionImpl.java Normal file
View File

@ -0,0 +1,435 @@
/*
* 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 org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathRuntimeException;
/**
* Calculates the Singular Value Decomposition of a matrix.
* <p>The Singular Value Decomposition of matrix A is a set of three matrices:
* U, &Sigma; and V such that A = U &times; &Sigma; &times; V<sup>T</sup>.
* Let A be an m &times; n matrix, then U is an m &times; m orthogonal matrix,
* &Sigma; is a m &times; n diagonal matrix with positive diagonal elements,
* and V is an n &times; n orthogonal matrix.</p>
*
* @version $Revision$ $Date$
* @since 2.0
*/
public class SingularValueDecompositionImpl implements SingularValueDecomposition {
/** Serializable version identifier. */
private static final long serialVersionUID = -2357152028714378552L;
/** Number of rows of the initial matrix. */
private int m;
/** Number of columns of the initial matrix. */
private int n;
/** Transformer to bidiagonal. */
private BiDiagonalTransformer transformer;
/** Main diagonal of the bidiagonal matrix. */
private double[] mainBidiagonal;
/** Secondary diagonal of the bidiagonal matrix. */
private double[] secondaryBidiagonal;
/** Main diagonal of the tridiagonal matrix. */
double[] mainTridiagonal;
/** Secondary diagonal of the tridiagonal matrix. */
double[] secondaryTridiagonal;
/** Eigen decomposition of the tridiagonal matrix. */
private EigenDecomposition eigenDecomposition;
/** Singular values. */
private double[] singularValues;
/** Cached value of U. */
private RealMatrix cachedU;
/** Cached value of U<sup>T</sup>. */
private RealMatrix cachedUt;
/** Cached value of S. */
private RealMatrix cachedS;
/** Cached value of V. */
private RealMatrix cachedV;
/** Cached value of V<sup>T</sup>. */
private RealMatrix cachedVt;
/**
* Build a new instance.
* <p>Note that {@link #decompose(RealMatrix)} <strong>must</strong> be called
* before any of the {@link #getU()}, {@link #getS()}, {@link #getV()},
* {@link #getSingularValues()}, {@link #getNorm()}, {@link #getConditionNumber()},
* {@link #getRank()}, {@link #solve(double[])}, {@link #solve(RealMatrix)},
* {@link #solve(RealVector)} or {@link #solve(RealVectorImpl)} methods can be
* called.</p>
* @see #decompose(RealMatrix)
*/
public SingularValueDecompositionImpl() {
}
/**
* Calculates the Singular Value Decomposition of the given matrix.
* <p>Calling this constructor is equivalent to first call the no-arguments
* constructor and then call {@link #decompose(RealMatrix)}.</p>
* @param matrix The matrix to decompose.
* @exception InvalidMatrixException (wrapping a {@link ConvergenceException}
* if algorithm fails to converge
*/
public SingularValueDecompositionImpl(RealMatrix matrix)
throws InvalidMatrixException {
decompose(matrix);
}
/** {@inheritDoc} */
public void decompose(final RealMatrix matrix) {
m = matrix.getRowDimension();
n = matrix.getColumnDimension();
cachedU = null;
cachedS = null;
cachedV = null;
cachedVt = null;
// transform the matrix to bidiagonal
transformer = new BiDiagonalTransformer(matrix);
mainBidiagonal = transformer.getMainDiagonalRef();
secondaryBidiagonal = transformer.getSecondaryDiagonalRef();
// compute Bt.B (if upper diagonal) or B.Bt (if lower diagonal)
mainTridiagonal = new double[mainBidiagonal.length];
secondaryTridiagonal = new double[mainBidiagonal.length - 1];
double a = mainBidiagonal[0];
mainTridiagonal[0] = a * a;
for (int i = 1; i < mainBidiagonal.length; ++i) {
final double b = secondaryBidiagonal[i - 1];
secondaryTridiagonal[i - 1] = a * b;
a = mainBidiagonal[i];
mainTridiagonal[i] = a * a + b * b;
}
// compute singular values
eigenDecomposition = new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal);
singularValues = eigenDecomposition.getEigenvalues();
for (int i = 0; i < singularValues.length; ++i) {
singularValues[i] = Math.sqrt(singularValues[i]);
}
}
/** {@inheritDoc} */
public RealMatrix getU()
throws InvalidMatrixException {
if (cachedU == null) {
checkDecomposed();
if (m >= n) {
// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
final double[][] eData = eigenDecomposition.getV().getData();
final double[][] iData = new double[m][];
double[] ei1 = eData[0];
iData[0] = ei1;
for (int i = 0; i < n - 1; ++i) {
// compute Bt.E.S^(-1) where E is the eigenvectors matrix
// we reuse the array from matrix E to store the result
final double[] ei0 = ei1;
ei1 = eData[i + 1];
iData[i + 1] = ei1;
for (int j = 0; j < n; ++j) {
ei0[j] = (mainBidiagonal[i] * ei0[j] +
secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
}
}
// last row
final double lastMain = mainBidiagonal[n - 1];
for (int j = 0; j < n; ++j) {
ei1[j] *= lastMain / singularValues[j];
}
for (int i = n; i < m; ++i) {
iData[i] = new double[n];
}
cachedU =
transformer.getU().multiply(new RealMatrixImpl(iData, false));
} else {
// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
cachedU = transformer.getU().multiply(eigenDecomposition.getV());
}
}
// return the cached matrix
return cachedU;
}
/** {@inheritDoc} */
public RealMatrix getUT()
throws InvalidMatrixException {
if (cachedUt == null) {
cachedUt = getU().transpose();
}
// return the cached matrix
return cachedUt;
}
/** {@inheritDoc} */
public RealMatrix getS()
throws InvalidMatrixException {
if (cachedS == null) {
checkDecomposed();
final int p = singularValues.length;
final double[][] sData = new double[p][p];
for (int i = 0; i < p; ++i) {
sData[i][i] = singularValues[i];
}
// cache the matrix for subsequent calls
cachedS = new RealMatrixImpl(sData, false);
}
return cachedS;
}
/** {@inheritDoc} */
public double[] getSingularValues()
throws InvalidMatrixException {
checkDecomposed();
return singularValues.clone();
}
/** {@inheritDoc} */
public RealMatrix getV()
throws InvalidMatrixException {
if (cachedV == null) {
checkDecomposed();
if (m >= n) {
// the tridiagonal matrix is Bt.B, where B is upper bidiagonal
cachedV = transformer.getV().multiply(eigenDecomposition.getV());
} else {
// the tridiagonal matrix is B.Bt, where B is lower bidiagonal
final double[][] eData = eigenDecomposition.getV().getData();
final double[][] iData = new double[n][];
double[] ei1 = eData[0];
iData[0] = ei1;
for (int i = 0; i < m - 1; ++i) {
// compute Bt.E.S^(-1) where E is the eigenvectors matrix
// we reuse the array from matrix E to store the result
final double[] ei0 = ei1;
ei1 = eData[i + 1];
iData[i + 1] = ei1;
for (int j = 0; j < m; ++j) {
ei0[j] = (mainBidiagonal[i] * ei0[j] +
secondaryBidiagonal[i] * ei1[j]) / singularValues[j];
}
}
// last row
final double lastMain = mainBidiagonal[m - 1];
for (int j = 0; j < m; ++j) {
ei1[j] *= lastMain / singularValues[j];
}
for (int i = m; i < n; ++i) {
iData[i] = new double[m];
}
cachedV =
transformer.getV().multiply(new RealMatrixImpl(iData, false));
}
}
// return the cached matrix
return cachedV;
}
/** {@inheritDoc} */
public RealMatrix getVT()
throws InvalidMatrixException {
if (cachedVt == null) {
cachedVt = getV().transpose();
}
// return the cached matrix
return cachedVt;
}
/** {@inheritDoc} */
public double getNorm()
throws InvalidMatrixException {
checkDecomposed();
return singularValues[0];
}
/** {@inheritDoc} */
public double getConditionNumber()
throws InvalidMatrixException {
checkDecomposed();
return singularValues[0] / singularValues[singularValues.length - 1];
}
/** {@inheritDoc} */
public int getRank()
throws IllegalStateException {
checkDecomposed();
final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]);
for (int i = singularValues.length - 1; i >= 0; --i) {
if (singularValues[i] > threshold) {
return i + 1;
}
}
return 0;
}
/** {@inheritDoc} */
public boolean isNonSingular()
throws IllegalStateException {
return getRank() == singularValues.length;
}
/** {@inheritDoc} */
public double[] solve(final double[] b)
throws IllegalArgumentException, InvalidMatrixException {
checkDecomposed();
if (b.length != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final double[] w = 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 getV().operate(w);
}
/** {@inheritDoc} */
public RealVector solve(final RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
try {
return solve((RealVectorImpl) b);
} catch (ClassCastException cce) {
checkDecomposed();
if (b.getDimension() != m) {
throw new IllegalArgumentException("constant vector has wrong length");
}
final RealVector w = 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 getV().operate(w);
}
}
/** Solve the linear equation A &times; X = B.
* <p>The A matrix is implicit here. It is </p>
* @param b right-hand side of the equation A &times; X = B
* @return a vector X such that A &times; X = B
* @throws IllegalArgumentException if matrices dimensions don't match
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public RealVectorImpl solve(final RealVectorImpl b)
throws IllegalArgumentException, InvalidMatrixException {
return new RealVectorImpl(solve(b.getDataRef()), false);
}
/** {@inheritDoc} */
public RealMatrix solve(final RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
checkDecomposed();
if (b.getRowDimension() != m) {
throw new IllegalArgumentException("Incorrect row dimension");
}
final RealMatrixImpl w = (RealMatrixImpl) 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 getV().multiply(w);
}
/** {@inheritDoc} */
public RealMatrix getInverse()
throws IllegalStateException, InvalidMatrixException {
checkDecomposed();
return solve(MatrixUtils.createRealIdentityMatrix(singularValues.length));
}
/**
* Check if {@link #decompose(RealMatrix)} has been called.
* @exception IllegalStateException if {@link #decompose(RealMatrix) decompose}
* has not been called
*/
private void checkDecomposed()
throws IllegalStateException {
if (singularValues == null) {
throw MathRuntimeException.createIllegalStateException("no matrix have been decomposed yet", null);
}
}
}

4
src/site/xdoc/changes.xml
View File

@ -78,8 +78,8 @@ The <action> type attribute can be add,update,fix,remove.
</action>
<action dev="luc" type="add" issue="MATH-220">
Added JAMA-like interfaces for eigen/singular problems. The implementation
of eigen-decomposition is based on the very quick dqd/dqds algorithms and some
parts of the MRRR algorithm. This leads to very fast and accurate solutions.
are based on the very quick dqd/dqds algorithms and some parts of the MRRR
algorithm. This leads to very fast and accurate solutions.
</action>
<action dev="luc" type="add" issue="MATH-220" >
Added JAMA-like interfaces for decomposition algorithms. These interfaces

332
src/test/org/apache/commons/math/linear/SingularValueDecompositionImplTest.java Executable file
View File

@ -0,0 +1,332 @@
/*
* 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;
public class SingularValueDecompositionImplTest extends TestCase {
private double[][] testSquare = {
{ 24.0 / 25.0, 43.0 / 25.0 },
{ 57.0 / 25.0, 24.0 / 25.0 }
};
private double[][] testNonSquare = {
{ -540.0 / 625.0, 963.0 / 625.0, -216.0 / 625.0 },
{ -1730.0 / 625.0, -744.0 / 625.0, 1008.0 / 625.0 },
{ -720.0 / 625.0, 1284.0 / 625.0, -288.0 / 625.0 },
{ -360.0 / 625.0, 192.0 / 625.0, 1756.0 / 625.0 },
};
private static final double normTolerance = 10e-14;
public SingularValueDecompositionImplTest(String name) {
super(name);
}
public static Test suite() {
TestSuite suite = new TestSuite(SingularValueDecompositionImplTest.class);
suite.setName("SingularValueDecompositionImpl Tests");
return suite;
}
public void testMoreRows() {
final double[] singularValues = { 123.456, 2.3, 1.001, 0.999 };
final int rows = singularValues.length + 2;
final int columns = singularValues.length;
Random r = new Random(15338437322523l);
SingularValueDecomposition svd =
new SingularValueDecompositionImpl(createTestMatrix(r, rows, columns, singularValues));
double[] computedSV = svd.getSingularValues();
assertEquals(singularValues.length, computedSV.length);
for (int i = 0; i < singularValues.length; ++i) {
assertEquals(singularValues[i], computedSV[i], 1.0e-10);
}
}
public void testMoreColumns() {
final double[] singularValues = { 123.456, 2.3, 1.001, 0.999 };
final int rows = singularValues.length;
final int columns = singularValues.length + 2;
Random r = new Random(732763225836210l);
SingularValueDecomposition svd =
new SingularValueDecompositionImpl(createTestMatrix(r, rows, columns, singularValues));
double[] computedSV = svd.getSingularValues();
assertEquals(singularValues.length, computedSV.length);
for (int i = 0; i < singularValues.length; ++i) {
assertEquals(singularValues[i], computedSV[i], 1.0e-10);
}
}
/** test dimensions */
public void testDimensions() {
RealMatrixImpl matrix = new RealMatrixImpl(testSquare, false);
final int m = matrix.getRowDimension();
final int n = matrix.getColumnDimension();
SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
assertEquals(m, svd.getU().getRowDimension());
assertEquals(m, svd.getU().getColumnDimension());
assertEquals(m, svd.getS().getColumnDimension());
assertEquals(n, svd.getS().getColumnDimension());
assertEquals(n, svd.getV().getRowDimension());
assertEquals(n, svd.getV().getColumnDimension());
}
/** test A = USVt */
public void testAEqualUSVt() {
checkAEqualUSVt(new RealMatrixImpl(testSquare, false));
checkAEqualUSVt(new RealMatrixImpl(testNonSquare, false));
checkAEqualUSVt(new RealMatrixImpl(testNonSquare, false).transpose());
}
public void checkAEqualUSVt(final RealMatrix matrix) {
SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
RealMatrix u = svd.getU();
RealMatrix s = svd.getS();
RealMatrix v = svd.getV();
double norm = u.multiply(s).multiply(v.transpose()).subtract(matrix).getNorm();
assertEquals(0, norm, normTolerance);
}
/** test that U is orthogonal */
public void testUOrthogonal() {
checkOrthogonal(new SingularValueDecompositionImpl(new RealMatrixImpl(testSquare, false)).getU());
checkOrthogonal(new SingularValueDecompositionImpl(new RealMatrixImpl(testNonSquare, false)).getU());
checkOrthogonal(new SingularValueDecompositionImpl(new RealMatrixImpl(testNonSquare, false).transpose()).getU());
}
/** test that V is orthogonal */
public void testVOrthogonal() {
checkOrthogonal(new SingularValueDecompositionImpl(new RealMatrixImpl(testSquare, false)).getV());
checkOrthogonal(new SingularValueDecompositionImpl(new RealMatrixImpl(testNonSquare, false)).getV());
checkOrthogonal(new SingularValueDecompositionImpl(new RealMatrixImpl(testNonSquare, false).transpose()).getV());
}
public void checkOrthogonal(final RealMatrix m) {
RealMatrix mTm = m.transpose().multiply(m);
RealMatrix id = MatrixUtils.createRealIdentityMatrix(mTm.getRowDimension());
assertEquals(0, mTm.subtract(id).getNorm(), normTolerance);
}
/** test solve dimension errors */
public void testSolveDimensionErrors() {
SingularValueDecomposition svd =
new SingularValueDecompositionImpl(new RealMatrixImpl(testSquare, false));
RealMatrix b = new RealMatrixImpl(new double[3][2]);
try {
svd.solve(b);
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
svd.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (IllegalArgumentException iae) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
svd.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() {
SingularValueDecomposition svd =
new SingularValueDecompositionImpl(new RealMatrixImpl(new double[][] {
{ 1.0, 0.0 },
{ 0.0, 0.0 }
}, false));
RealMatrix b = new RealMatrixImpl(new double[2][2]);
try {
svd.solve(b);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
svd.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
svd.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
// expected behavior
} catch (Exception e) {
fail("wrong exception caught");
}
try {
svd.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() {
SingularValueDecomposition svd =
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, svd.solve(b).subtract(xRef).getNorm(), normTolerance);
// using double[]
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
new RealVectorImpl(svd.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
// using RealMatrixImpl
for (int i = 0; i < b.getColumnDimension(); ++i) {
assertEquals(0,
svd.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,
svd.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
1.0e-13);
}
}
/** test matrices values */
public void testMatricesValues1() {
SingularValueDecomposition svd =
new SingularValueDecompositionImpl(new RealMatrixImpl(testSquare, false));
RealMatrix uRef = new RealMatrixImpl(new double[][] {
{ 3.0 / 5.0, -4.0 / 5.0 },
{ 4.0 / 5.0, 3.0 / 5.0 }
});
RealMatrix sRef = new RealMatrixImpl(new double[][] {
{ 3.0, 0.0 },
{ 0.0, 1.0 }
});
RealMatrix vRef = new RealMatrixImpl(new double[][] {
{ 4.0 / 5.0, 3.0 / 5.0 },
{ 3.0 / 5.0, -4.0 / 5.0 }
});
// check values against known references
RealMatrix u = svd.getU();
assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
RealMatrix s = svd.getS();
assertEquals(0, s.subtract(sRef).getNorm(), normTolerance);
RealMatrix v = svd.getV();
assertEquals(0, v.subtract(vRef).getNorm(), normTolerance);
// check the same cached instance is returned the second time
assertTrue(u == svd.getU());
assertTrue(s == svd.getS());
assertTrue(v == svd.getV());
}
/** test matrices values */
public void testMatricesValues2() {
RealMatrix uRef = new RealMatrixImpl(new double[][] {
{ 0.0 / 5.0, 3.0 / 5.0, 0.0 / 5.0 },
{ -4.0 / 5.0, 0.0 / 5.0, -3.0 / 5.0 },
{ 0.0 / 5.0, 4.0 / 5.0, 0.0 / 5.0 },
{ -3.0 / 5.0, 0.0 / 5.0, 4.0 / 5.0 }
});
RealMatrix sRef = new RealMatrixImpl(new double[][] {
{ 4.0, 0.0, 0.0 },
{ 0.0, 3.0, 0.0 },
{ 0.0, 0.0, 2.0 }
});
RealMatrix vRef = new RealMatrixImpl(new double[][] {
{ 80.0 / 125.0, -60.0 / 125.0, 75.0 / 125.0 },
{ 24.0 / 125.0, 107.0 / 125.0, 60.0 / 125.0 },
{ -93.0 / 125.0, -24.0 / 125.0, 80.0 / 125.0 }
});
// check values against known references
SingularValueDecomposition svd = new SingularValueDecompositionImpl();
svd.decompose(new RealMatrixImpl(testNonSquare, false));
RealMatrix u = svd.getU();
assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
RealMatrix s = svd.getS();
assertEquals(0, s.subtract(sRef).getNorm(), normTolerance);
RealMatrix v = svd.getV();
assertEquals(0, v.subtract(vRef).getNorm(), normTolerance);
// check the same cached instance is returned the second time
assertTrue(u == svd.getU());
assertTrue(s == svd.getS());
assertTrue(v == svd.getV());
}
/** test condition number */
public void testConditionNumber() {
SingularValueDecompositionImpl svd =
new SingularValueDecompositionImpl(new RealMatrixImpl(testSquare, false));
assertEquals(3.0, svd.getConditionNumber(), 1.0e-15);
}
private RealMatrix createTestMatrix(final Random r, final int rows, final int columns,
final double[] singularValues) {
final RealMatrix u =
EigenDecompositionImplTest.createOrthogonalMatrix(r, rows);
final RealMatrix d =
EigenDecompositionImplTest.createDiagonalMatrix(singularValues, rows, columns);
final RealMatrix v =
EigenDecompositionImplTest.createOrthogonalMatrix(r, columns);
return u.multiply(d).multiply(v);
}
}