Merged SingularValueDecomposition and SingularValueDecompositionImpl (see MATH-662).
git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1175108 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Sebastien Brisard
parent
b4764661a3
commit
39d86e32e2
|
|
@ -14,14 +14,15 @@
|
|||
* 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.exception.NumberIsTooLargeException;
|
||||
import org.apache.commons.math.exception.util.LocalizedFormats;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
import org.apache.commons.math.util.MathUtils;
|
||||
|
||||
/**
|
||||
* An interface to classes that implement an algorithm to calculate the
|
||||
* Singular Value Decomposition of a real matrix.
|
||||
* Calculates the compact Singular Value Decomposition of a matrix.
|
||||
* <p>
|
||||
* The Singular Value Decomposition of matrix A is a set of three matrices: U,
|
||||
* Σ and V such that A = U × Σ × V<sup>T</sup>. Let A be
|
||||
|
|
@ -30,15 +31,15 @@ package org.apache.commons.math.linear;
|
|||
* n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where
|
||||
* p=min(m,n).
|
||||
* </p>
|
||||
* <p>This interface is similar to the class with similar name from the
|
||||
* <p>This class is similar to the class with similar name from the
|
||||
* <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library, with the
|
||||
* following changes:</p>
|
||||
* <ul>
|
||||
* <li>the <code>norm2</code> method which has been renamed as {@link #getNorm()
|
||||
* <li>the {@code norm2} method which has been renamed as {@link #getNorm()
|
||||
* getNorm},</li>
|
||||
* <li>the <code>cond</code> method which has been renamed as {@link
|
||||
* <li>the {@code cond} method which has been renamed as {@link
|
||||
* #getConditionNumber() getConditionNumber},</li>
|
||||
* <li>the <code>rank</code> method which has been renamed as {@link #getRank()
|
||||
* <li>the {@code rank} method which has been renamed as {@link #getRank()
|
||||
* getRank},</li>
|
||||
* <li>a {@link #getUT() getUT} method has been added,</li>
|
||||
* <li>a {@link #getVT() getVT} method has been added,</li>
|
||||
|
|
@ -48,9 +49,432 @@ package org.apache.commons.math.linear;
|
|||
* @see <a href="http://mathworld.wolfram.com/SingularValueDecomposition.html">MathWorld</a>
|
||||
* @see <a href="http://en.wikipedia.org/wiki/Singular_value_decomposition">Wikipedia</a>
|
||||
* @version $Id$
|
||||
* @since 2.0
|
||||
* @since 2.0 (changed to concrete class in 3.0)
|
||||
*/
|
||||
public interface SingularValueDecomposition {
|
||||
public class SingularValueDecomposition {
|
||||
/** Relative threshold for small singular values. */
|
||||
private static final double EPS = 0x1.0p-52;
|
||||
/** Absolute threshold for small singular values. */
|
||||
private static final double TINY = 0x1.0p-966;
|
||||
/** Computed singular values. */
|
||||
private final double[] singularValues;
|
||||
/** max(row dimension, column dimension). */
|
||||
private final int m;
|
||||
/** min(row dimension, column dimension). */
|
||||
private final int n;
|
||||
/** Indicator for transposed matrix. */
|
||||
private final boolean transposed;
|
||||
/** Cached value of U matrix. */
|
||||
private final RealMatrix cachedU;
|
||||
/** Cached value of transposed U matrix. */
|
||||
private RealMatrix cachedUt;
|
||||
/** Cached value of S (diagonal) matrix. */
|
||||
private RealMatrix cachedS;
|
||||
/** Cached value of V matrix. */
|
||||
private final RealMatrix cachedV;
|
||||
/** Cached value of transposed V matrix. */
|
||||
private RealMatrix cachedVt;
|
||||
/**
|
||||
* Tolerance value for small singular values, calculated once we have
|
||||
* populated "singularValues".
|
||||
**/
|
||||
private final double tol;
|
||||
|
||||
/**
|
||||
* Calculates the compact Singular Value Decomposition of the given matrix.
|
||||
*
|
||||
* @param matrix Matrix to decompose.
|
||||
*/
|
||||
public SingularValueDecomposition(final RealMatrix matrix) {
|
||||
final double[][] A;
|
||||
|
||||
// "m" is always the largest dimension.
|
||||
if (matrix.getRowDimension() < matrix.getColumnDimension()) {
|
||||
transposed = true;
|
||||
A = matrix.transpose().getData();
|
||||
m = matrix.getColumnDimension();
|
||||
n = matrix.getRowDimension();
|
||||
} else {
|
||||
transposed = false;
|
||||
A = matrix.getData();
|
||||
m = matrix.getRowDimension();
|
||||
n = matrix.getColumnDimension();
|
||||
}
|
||||
|
||||
singularValues = new double[n];
|
||||
final double[][] U = new double[m][n];
|
||||
final double[][] V = new double[n][n];
|
||||
final double[] e = new double[n];
|
||||
final double[] work = new double[m];
|
||||
// Reduce A to bidiagonal form, storing the diagonal elements
|
||||
// in s and the super-diagonal elements in e.
|
||||
final int nct = FastMath.min(m - 1, n);
|
||||
final int nrt = FastMath.max(0, n - 2);
|
||||
for (int k = 0; k < FastMath.max(nct, nrt); k++) {
|
||||
if (k < nct) {
|
||||
// Compute the transformation for the k-th column and
|
||||
// place the k-th diagonal in s[k].
|
||||
// Compute 2-norm of k-th column without under/overflow.
|
||||
singularValues[k] = 0;
|
||||
for (int i = k; i < m; i++) {
|
||||
singularValues[k] = FastMath.hypot(singularValues[k], A[i][k]);
|
||||
}
|
||||
if (singularValues[k] != 0) {
|
||||
if (A[k][k] < 0) {
|
||||
singularValues[k] = -singularValues[k];
|
||||
}
|
||||
for (int i = k; i < m; i++) {
|
||||
A[i][k] /= singularValues[k];
|
||||
}
|
||||
A[k][k] += 1;
|
||||
}
|
||||
singularValues[k] = -singularValues[k];
|
||||
}
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
if (k < nct &&
|
||||
singularValues[k] != 0) {
|
||||
// Apply the transformation.
|
||||
double t = 0;
|
||||
for (int i = k; i < m; i++) {
|
||||
t += A[i][k] * A[i][j];
|
||||
}
|
||||
t = -t / A[k][k];
|
||||
for (int i = k; i < m; i++) {
|
||||
A[i][j] += t * A[i][k];
|
||||
}
|
||||
}
|
||||
// Place the k-th row of A into e for the
|
||||
// subsequent calculation of the row transformation.
|
||||
e[j] = A[k][j];
|
||||
}
|
||||
if (k < nct) {
|
||||
// Place the transformation in U for subsequent back
|
||||
// multiplication.
|
||||
for (int i = k; i < m; i++) {
|
||||
U[i][k] = A[i][k];
|
||||
}
|
||||
}
|
||||
if (k < nrt) {
|
||||
// Compute the k-th row transformation and place the
|
||||
// k-th super-diagonal in e[k].
|
||||
// Compute 2-norm without under/overflow.
|
||||
e[k] = 0;
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
e[k] = FastMath.hypot(e[k], e[i]);
|
||||
}
|
||||
if (e[k] != 0) {
|
||||
if (e[k + 1] < 0) {
|
||||
e[k] = -e[k];
|
||||
}
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
e[i] /= e[k];
|
||||
}
|
||||
e[k + 1] += 1;
|
||||
}
|
||||
e[k] = -e[k];
|
||||
if (k + 1 < m &&
|
||||
e[k] != 0) {
|
||||
// Apply the transformation.
|
||||
for (int i = k + 1; i < m; i++) {
|
||||
work[i] = 0;
|
||||
}
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
for (int i = k + 1; i < m; i++) {
|
||||
work[i] += e[j] * A[i][j];
|
||||
}
|
||||
}
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
final double t = -e[j] / e[k + 1];
|
||||
for (int i = k + 1; i < m; i++) {
|
||||
A[i][j] += t * work[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Place the transformation in V for subsequent
|
||||
// back multiplication.
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
V[i][k] = e[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
// Set up the final bidiagonal matrix or order p.
|
||||
int p = n;
|
||||
if (nct < n) {
|
||||
singularValues[nct] = A[nct][nct];
|
||||
}
|
||||
if (m < p) {
|
||||
singularValues[p - 1] = 0;
|
||||
}
|
||||
if (nrt + 1 < p) {
|
||||
e[nrt] = A[nrt][p - 1];
|
||||
}
|
||||
e[p - 1] = 0;
|
||||
|
||||
// Generate U.
|
||||
for (int j = nct; j < n; j++) {
|
||||
for (int i = 0; i < m; i++) {
|
||||
U[i][j] = 0;
|
||||
}
|
||||
U[j][j] = 1;
|
||||
}
|
||||
for (int k = nct - 1; k >= 0; k--) {
|
||||
if (singularValues[k] != 0) {
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
double t = 0;
|
||||
for (int i = k; i < m; i++) {
|
||||
t += U[i][k] * U[i][j];
|
||||
}
|
||||
t = -t / U[k][k];
|
||||
for (int i = k; i < m; i++) {
|
||||
U[i][j] += t * U[i][k];
|
||||
}
|
||||
}
|
||||
for (int i = k; i < m; i++) {
|
||||
U[i][k] = -U[i][k];
|
||||
}
|
||||
U[k][k] = 1 + U[k][k];
|
||||
for (int i = 0; i < k - 1; i++) {
|
||||
U[i][k] = 0;
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < m; i++) {
|
||||
U[i][k] = 0;
|
||||
}
|
||||
U[k][k] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// Generate V.
|
||||
for (int k = n - 1; k >= 0; k--) {
|
||||
if (k < nrt &&
|
||||
e[k] != 0) {
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
double t = 0;
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
t += V[i][k] * V[i][j];
|
||||
}
|
||||
t = -t / V[k + 1][k];
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
V[i][j] += t * V[i][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < n; i++) {
|
||||
V[i][k] = 0;
|
||||
}
|
||||
V[k][k] = 1;
|
||||
}
|
||||
|
||||
// Main iteration loop for the singular values.
|
||||
final int pp = p - 1;
|
||||
int iter = 0;
|
||||
while (p > 0) {
|
||||
int k;
|
||||
int kase;
|
||||
// Here is where a test for too many iterations would go.
|
||||
// This section of the program inspects for
|
||||
// negligible elements in the s and e arrays. On
|
||||
// completion the variables kase and k are set as follows.
|
||||
// kase = 1 if s(p) and e[k-1] are negligible and k<p
|
||||
// kase = 2 if s(k) is negligible and k<p
|
||||
// kase = 3 if e[k-1] is negligible, k<p, and
|
||||
// s(k), ..., s(p) are not negligible (qr step).
|
||||
// kase = 4 if e(p-1) is negligible (convergence).
|
||||
for (k = p - 2; k >= 0; k--) {
|
||||
final double threshold
|
||||
= TINY + EPS * (FastMath.abs(singularValues[k]) +
|
||||
FastMath.abs(singularValues[k + 1]));
|
||||
if (FastMath.abs(e[k]) <= threshold) {
|
||||
e[k] = 0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (k == p - 2) {
|
||||
kase = 4;
|
||||
} else {
|
||||
int ks;
|
||||
for (ks = p - 1; ks >= k; ks--) {
|
||||
if (ks == k) {
|
||||
break;
|
||||
}
|
||||
final double t = (ks != p ? FastMath.abs(e[ks]) : 0) +
|
||||
(ks != k + 1 ? FastMath.abs(e[ks - 1]) : 0);
|
||||
if (FastMath.abs(singularValues[ks]) <= TINY + EPS * t) {
|
||||
singularValues[ks] = 0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (ks == k) {
|
||||
kase = 3;
|
||||
} else if (ks == p - 1) {
|
||||
kase = 1;
|
||||
} else {
|
||||
kase = 2;
|
||||
k = ks;
|
||||
}
|
||||
}
|
||||
k++;
|
||||
// Perform the task indicated by kase.
|
||||
switch (kase) {
|
||||
// Deflate negligible s(p).
|
||||
case 1: {
|
||||
double f = e[p - 2];
|
||||
e[p - 2] = 0;
|
||||
for (int j = p - 2; j >= k; j--) {
|
||||
double t = FastMath.hypot(singularValues[j], f);
|
||||
final double cs = singularValues[j] / t;
|
||||
final double sn = f / t;
|
||||
singularValues[j] = t;
|
||||
if (j != k) {
|
||||
f = -sn * e[j - 1];
|
||||
e[j - 1] = cs * e[j - 1];
|
||||
}
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
t = cs * V[i][j] + sn * V[i][p - 1];
|
||||
V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1];
|
||||
V[i][j] = t;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
// Split at negligible s(k).
|
||||
case 2: {
|
||||
double f = e[k - 1];
|
||||
e[k - 1] = 0;
|
||||
for (int j = k; j < p; j++) {
|
||||
double t = FastMath.hypot(singularValues[j], f);
|
||||
final double cs = singularValues[j] / t;
|
||||
final double sn = f / t;
|
||||
singularValues[j] = t;
|
||||
f = -sn * e[j];
|
||||
e[j] = cs * e[j];
|
||||
|
||||
for (int i = 0; i < m; i++) {
|
||||
t = cs * U[i][j] + sn * U[i][k - 1];
|
||||
U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1];
|
||||
U[i][j] = t;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
// Perform one qr step.
|
||||
case 3: {
|
||||
// Calculate the shift.
|
||||
final double scale = FastMath.max(FastMath.max(FastMath.max(FastMath.max(
|
||||
FastMath.abs(singularValues[p - 1]), FastMath.abs(singularValues[p - 2])), FastMath.abs(e[p - 2])),
|
||||
FastMath.abs(singularValues[k])), FastMath.abs(e[k]));
|
||||
final double sp = singularValues[p - 1] / scale;
|
||||
final double spm1 = singularValues[p - 2] / scale;
|
||||
final double epm1 = e[p - 2] / scale;
|
||||
final double sk = singularValues[k] / scale;
|
||||
final double ek = e[k] / scale;
|
||||
final double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
|
||||
final double c = (sp * epm1) * (sp * epm1);
|
||||
double shift = 0;
|
||||
if (b != 0 ||
|
||||
c != 0) {
|
||||
shift = FastMath.sqrt(b * b + c);
|
||||
if (b < 0) {
|
||||
shift = -shift;
|
||||
}
|
||||
shift = c / (b + shift);
|
||||
}
|
||||
double f = (sk + sp) * (sk - sp) + shift;
|
||||
double g = sk * ek;
|
||||
// Chase zeros.
|
||||
for (int j = k; j < p - 1; j++) {
|
||||
double t = FastMath.hypot(f, g);
|
||||
double cs = f / t;
|
||||
double sn = g / t;
|
||||
if (j != k) {
|
||||
e[j - 1] = t;
|
||||
}
|
||||
f = cs * singularValues[j] + sn * e[j];
|
||||
e[j] = cs * e[j] - sn * singularValues[j];
|
||||
g = sn * singularValues[j + 1];
|
||||
singularValues[j + 1] = cs * singularValues[j + 1];
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
t = cs * V[i][j] + sn * V[i][j + 1];
|
||||
V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1];
|
||||
V[i][j] = t;
|
||||
}
|
||||
t = FastMath.hypot(f, g);
|
||||
cs = f / t;
|
||||
sn = g / t;
|
||||
singularValues[j] = t;
|
||||
f = cs * e[j] + sn * singularValues[j + 1];
|
||||
singularValues[j + 1] = -sn * e[j] + cs * singularValues[j + 1];
|
||||
g = sn * e[j + 1];
|
||||
e[j + 1] = cs * e[j + 1];
|
||||
if (j < m - 1) {
|
||||
for (int i = 0; i < m; i++) {
|
||||
t = cs * U[i][j] + sn * U[i][j + 1];
|
||||
U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1];
|
||||
U[i][j] = t;
|
||||
}
|
||||
}
|
||||
}
|
||||
e[p - 2] = f;
|
||||
iter = iter + 1;
|
||||
}
|
||||
break;
|
||||
// Convergence.
|
||||
default: {
|
||||
// Make the singular values positive.
|
||||
if (singularValues[k] <= 0) {
|
||||
singularValues[k] = singularValues[k] < 0 ? -singularValues[k] : 0;
|
||||
|
||||
for (int i = 0; i <= pp; i++) {
|
||||
V[i][k] = -V[i][k];
|
||||
}
|
||||
}
|
||||
// Order the singular values.
|
||||
while (k < pp) {
|
||||
if (singularValues[k] >= singularValues[k + 1]) {
|
||||
break;
|
||||
}
|
||||
double t = singularValues[k];
|
||||
singularValues[k] = singularValues[k + 1];
|
||||
singularValues[k + 1] = t;
|
||||
if (k < n - 1) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
t = V[i][k + 1];
|
||||
V[i][k + 1] = V[i][k];
|
||||
V[i][k] = t;
|
||||
}
|
||||
}
|
||||
if (k < m - 1) {
|
||||
for (int i = 0; i < m; i++) {
|
||||
t = U[i][k + 1];
|
||||
U[i][k + 1] = U[i][k];
|
||||
U[i][k] = t;
|
||||
}
|
||||
}
|
||||
k++;
|
||||
}
|
||||
iter = 0;
|
||||
p--;
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// Set the small value tolerance used to calculate rank and pseudo-inverse
|
||||
tol = FastMath.max(m * singularValues[0] * EPS,
|
||||
FastMath.sqrt(MathUtils.SAFE_MIN));
|
||||
|
||||
if (!transposed) {
|
||||
cachedU = MatrixUtils.createRealMatrix(U);
|
||||
cachedV = MatrixUtils.createRealMatrix(V);
|
||||
} else {
|
||||
cachedU = MatrixUtils.createRealMatrix(V);
|
||||
cachedV = MatrixUtils.createRealMatrix(U);
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the matrix U of the decomposition.
|
||||
|
|
@ -58,7 +482,11 @@ public interface SingularValueDecomposition {
|
|||
* @return the U matrix
|
||||
* @see #getUT()
|
||||
*/
|
||||
RealMatrix getU();
|
||||
public RealMatrix getU() {
|
||||
// return the cached matrix
|
||||
return cachedU;
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the transpose of the matrix U of the decomposition.
|
||||
|
|
@ -66,7 +494,13 @@ public interface SingularValueDecomposition {
|
|||
* @return the U matrix (or null if decomposed matrix is singular)
|
||||
* @see #getU()
|
||||
*/
|
||||
RealMatrix getUT();
|
||||
public RealMatrix getUT() {
|
||||
if (cachedUt == null) {
|
||||
cachedUt = getU().transpose();
|
||||
}
|
||||
// return the cached matrix
|
||||
return cachedUt;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the diagonal matrix Σ of the decomposition.
|
||||
|
|
@ -74,7 +508,13 @@ public interface SingularValueDecomposition {
|
|||
* non-increasing order, for compatibility with Jama.</p>
|
||||
* @return the Σ matrix
|
||||
*/
|
||||
RealMatrix getS();
|
||||
public RealMatrix getS() {
|
||||
if (cachedS == null) {
|
||||
// cache the matrix for subsequent calls
|
||||
cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
|
||||
}
|
||||
return cachedS;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the diagonal elements of the matrix Σ of the decomposition.
|
||||
|
|
@ -82,7 +522,9 @@ public interface SingularValueDecomposition {
|
|||
* compatibility with Jama.</p>
|
||||
* @return the diagonal elements of the Σ matrix
|
||||
*/
|
||||
double[] getSingularValues();
|
||||
public double[] getSingularValues() {
|
||||
return singularValues.clone();
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the matrix V of the decomposition.
|
||||
|
|
@ -90,7 +532,10 @@ public interface SingularValueDecomposition {
|
|||
* @return the V matrix (or null if decomposed matrix is singular)
|
||||
* @see #getVT()
|
||||
*/
|
||||
RealMatrix getV();
|
||||
public RealMatrix getV() {
|
||||
// return the cached matrix
|
||||
return cachedV;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the transpose of the matrix V of the decomposition.
|
||||
|
|
@ -98,7 +543,13 @@ public interface SingularValueDecomposition {
|
|||
* @return the V matrix (or null if decomposed matrix is singular)
|
||||
* @see #getV()
|
||||
*/
|
||||
RealMatrix getVT();
|
||||
public RealMatrix getVT() {
|
||||
if (cachedVt == null) {
|
||||
cachedVt = getV().transpose();
|
||||
}
|
||||
// return the cached matrix
|
||||
return cachedVt;
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the n × n covariance matrix.
|
||||
|
|
@ -111,7 +562,33 @@ public interface SingularValueDecomposition {
|
|||
* @exception IllegalArgumentException if minSingularValue is larger than
|
||||
* the largest singular value, meaning all singular values are ignored
|
||||
*/
|
||||
RealMatrix getCovariance(double minSingularValue);
|
||||
public RealMatrix getCovariance(final double minSingularValue) {
|
||||
// get the number of singular values to consider
|
||||
final int p = singularValues.length;
|
||||
int dimension = 0;
|
||||
while (dimension < p &&
|
||||
singularValues[dimension] >= minSingularValue) {
|
||||
++dimension;
|
||||
}
|
||||
|
||||
if (dimension == 0) {
|
||||
throw new NumberIsTooLargeException(LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE,
|
||||
minSingularValue, singularValues[0], true);
|
||||
}
|
||||
|
||||
final double[][] data = new double[dimension][p];
|
||||
getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
public void visit(final int row, final int column,
|
||||
final double value) {
|
||||
data[row][column] = value / singularValues[row];
|
||||
}
|
||||
}, 0, dimension - 1, 0, p - 1);
|
||||
|
||||
RealMatrix jv = new Array2DRowRealMatrix(data, false);
|
||||
return jv.transpose().multiply(jv);
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the L<sub>2</sub> norm of the matrix.
|
||||
|
|
@ -120,13 +597,28 @@ public interface SingularValueDecomposition {
|
|||
* (i.e. the traditional euclidian norm).</p>
|
||||
* @return norm
|
||||
*/
|
||||
double getNorm();
|
||||
public double getNorm() {
|
||||
return singularValues[0];
|
||||
}
|
||||
|
||||
/**
|
||||
* Return the condition number of the matrix.
|
||||
* @return condition number of the matrix
|
||||
*/
|
||||
double getConditionNumber();
|
||||
public double getConditionNumber() {
|
||||
return singularValues[0] / singularValues[n - 1];
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes the inverse of the condition number.
|
||||
* In cases of rank deficiency, the {@link #getConditionNumber() condition
|
||||
* number} will become undefined.
|
||||
*
|
||||
* @return the inverse of the condition number.
|
||||
*/
|
||||
public double getInverseConditionNumber() {
|
||||
return singularValues[n - 1] / singularValues[0];
|
||||
}
|
||||
|
||||
/**
|
||||
* Return the effective numerical matrix rank.
|
||||
|
|
@ -136,12 +628,106 @@ public interface SingularValueDecomposition {
|
|||
* is the least significant bit of the largest singular value.</p>
|
||||
* @return effective numerical matrix rank
|
||||
*/
|
||||
int getRank();
|
||||
public int getRank() {
|
||||
int r = 0;
|
||||
for (int i = 0; i < singularValues.length; i++) {
|
||||
if (singularValues[i] > tol) {
|
||||
r++;
|
||||
}
|
||||
}
|
||||
return r;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get a solver for finding the A × X = B solution in least square sense.
|
||||
* @return a solver
|
||||
*/
|
||||
DecompositionSolver getSolver();
|
||||
public DecompositionSolver getSolver() {
|
||||
return new Solver(singularValues, getUT(), getV(), getRank() == m, tol);
|
||||
}
|
||||
|
||||
/** Specialized solver. */
|
||||
private static class Solver implements DecompositionSolver {
|
||||
/** Pseudo-inverse of the initial matrix. */
|
||||
private final RealMatrix pseudoInverse;
|
||||
/** Singularity indicator. */
|
||||
private boolean nonSingular;
|
||||
|
||||
/**
|
||||
* Build a solver from decomposed matrix.
|
||||
*
|
||||
* @param singularValues Singular values.
|
||||
* @param uT U<sup>T</sup> matrix of the decomposition.
|
||||
* @param v V matrix of the decomposition.
|
||||
* @param nonSingular Singularity indicator.
|
||||
* @param tol tolerance for singular values
|
||||
*/
|
||||
private Solver(final double[] singularValues, final RealMatrix uT,
|
||||
final RealMatrix v, final boolean nonSingular, final double tol) {
|
||||
final double[][] suT = uT.getData();
|
||||
for (int i = 0; i < singularValues.length; ++i) {
|
||||
final double a;
|
||||
if (singularValues[i] > tol) {
|
||||
a = 1 / singularValues[i];
|
||||
} else {
|
||||
a = 0;
|
||||
}
|
||||
final double[] suTi = suT[i];
|
||||
for (int j = 0; j < suTi.length; ++j) {
|
||||
suTi[j] *= a;
|
||||
}
|
||||
}
|
||||
pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false));
|
||||
this.nonSingular = nonSingular;
|
||||
}
|
||||
|
||||
/**
|
||||
* Solve the linear equation A × X = B in least square sense.
|
||||
* <p>
|
||||
* The m×n matrix A may not be square, the solution X is such that
|
||||
* ||A × X - B|| is minimal.
|
||||
* </p>
|
||||
* @param b Right-hand side of the equation A × X = B
|
||||
* @return a vector X that minimizes the two norm of A × X - B
|
||||
* @throws org.apache.commons.math.exception.DimensionMismatchException
|
||||
* if the matrices dimensions do not match.
|
||||
*/
|
||||
public RealVector solve(final RealVector b) {
|
||||
return pseudoInverse.operate(b);
|
||||
}
|
||||
|
||||
/**
|
||||
* Solve the linear equation A × X = B in least square sense.
|
||||
* <p>
|
||||
* The m×n matrix A may not be square, the solution X is such that
|
||||
* ||A × X - B|| is minimal.
|
||||
* </p>
|
||||
*
|
||||
* @param b Right-hand side of the equation A × X = B
|
||||
* @return a matrix X that minimizes the two norm of A × X - B
|
||||
* @throws org.apache.commons.math.exception.DimensionMismatchException
|
||||
* if the matrices dimensions do not match.
|
||||
*/
|
||||
public RealMatrix solve(final RealMatrix b) {
|
||||
return pseudoInverse.multiply(b);
|
||||
}
|
||||
|
||||
/**
|
||||
* Check if the decomposed matrix is non-singular.
|
||||
*
|
||||
* @return {@code true} if the decomposed matrix is non-singular.
|
||||
*/
|
||||
public boolean isNonSingular() {
|
||||
return nonSingular;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get the pseudo-inverse of the decomposed matrix.
|
||||
*
|
||||
* @return the inverse matrix.
|
||||
*/
|
||||
public RealMatrix getInverse() {
|
||||
return pseudoInverse;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,657 +0,0 @@
|
|||
/*
|
||||
* 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.exception.NumberIsTooLargeException;
|
||||
import org.apache.commons.math.exception.util.LocalizedFormats;
|
||||
import org.apache.commons.math.util.FastMath;
|
||||
import org.apache.commons.math.util.MathUtils;
|
||||
|
||||
/**
|
||||
* Calculates the compact Singular Value Decomposition of a matrix.
|
||||
* <p>
|
||||
* The Singular Value Decomposition of matrix A is a set of three matrices: U,
|
||||
* Σ and V such that A = U × Σ × V<sup>T</sup>. Let A be
|
||||
* a m × n matrix, then U is a m × p orthogonal matrix, Σ is a
|
||||
* p × p diagonal matrix with positive or null elements, V is a p ×
|
||||
* n orthogonal matrix (hence V<sup>T</sup> is also orthogonal) where
|
||||
* p=min(m,n).
|
||||
* </p>
|
||||
* @version $Id$
|
||||
* @since 2.0
|
||||
*/
|
||||
public class SingularValueDecompositionImpl implements SingularValueDecomposition {
|
||||
/** Relative threshold for small singular values. */
|
||||
private static final double EPS = 0x1.0p-52;
|
||||
/** Absolute threshold for small singular values. */
|
||||
private static final double TINY = 0x1.0p-966;
|
||||
/** Computed singular values. */
|
||||
private final double[] singularValues;
|
||||
/** max(row dimension, column dimension). */
|
||||
private final int m;
|
||||
/** min(row dimension, column dimension). */
|
||||
private final int n;
|
||||
/** Indicator for transposed matrix. */
|
||||
private final boolean transposed;
|
||||
/** Cached value of U matrix. */
|
||||
private final RealMatrix cachedU;
|
||||
/** Cached value of transposed U matrix. */
|
||||
private RealMatrix cachedUt;
|
||||
/** Cached value of S (diagonal) matrix. */
|
||||
private RealMatrix cachedS;
|
||||
/** Cached value of V matrix. */
|
||||
private final RealMatrix cachedV;
|
||||
/** Cached value of transposed V matrix. */
|
||||
private RealMatrix cachedVt;
|
||||
/**
|
||||
* Tolerance value for small singular values, calculated once we have
|
||||
* populated "singularValues".
|
||||
**/
|
||||
private final double tol;
|
||||
|
||||
/**
|
||||
* Calculates the compact Singular Value Decomposition of the given matrix.
|
||||
*
|
||||
* @param matrix Matrix to decompose.
|
||||
*/
|
||||
public SingularValueDecompositionImpl(final RealMatrix matrix) {
|
||||
final double[][] A;
|
||||
|
||||
// "m" is always the largest dimension.
|
||||
if (matrix.getRowDimension() < matrix.getColumnDimension()) {
|
||||
transposed = true;
|
||||
A = matrix.transpose().getData();
|
||||
m = matrix.getColumnDimension();
|
||||
n = matrix.getRowDimension();
|
||||
} else {
|
||||
transposed = false;
|
||||
A = matrix.getData();
|
||||
m = matrix.getRowDimension();
|
||||
n = matrix.getColumnDimension();
|
||||
}
|
||||
|
||||
singularValues = new double[n];
|
||||
final double[][] U = new double[m][n];
|
||||
final double[][] V = new double[n][n];
|
||||
final double[] e = new double[n];
|
||||
final double[] work = new double[m];
|
||||
// Reduce A to bidiagonal form, storing the diagonal elements
|
||||
// in s and the super-diagonal elements in e.
|
||||
final int nct = FastMath.min(m - 1, n);
|
||||
final int nrt = FastMath.max(0, n - 2);
|
||||
for (int k = 0; k < FastMath.max(nct, nrt); k++) {
|
||||
if (k < nct) {
|
||||
// Compute the transformation for the k-th column and
|
||||
// place the k-th diagonal in s[k].
|
||||
// Compute 2-norm of k-th column without under/overflow.
|
||||
singularValues[k] = 0;
|
||||
for (int i = k; i < m; i++) {
|
||||
singularValues[k] = FastMath.hypot(singularValues[k], A[i][k]);
|
||||
}
|
||||
if (singularValues[k] != 0) {
|
||||
if (A[k][k] < 0) {
|
||||
singularValues[k] = -singularValues[k];
|
||||
}
|
||||
for (int i = k; i < m; i++) {
|
||||
A[i][k] /= singularValues[k];
|
||||
}
|
||||
A[k][k] += 1;
|
||||
}
|
||||
singularValues[k] = -singularValues[k];
|
||||
}
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
if (k < nct &&
|
||||
singularValues[k] != 0) {
|
||||
// Apply the transformation.
|
||||
double t = 0;
|
||||
for (int i = k; i < m; i++) {
|
||||
t += A[i][k] * A[i][j];
|
||||
}
|
||||
t = -t / A[k][k];
|
||||
for (int i = k; i < m; i++) {
|
||||
A[i][j] += t * A[i][k];
|
||||
}
|
||||
}
|
||||
// Place the k-th row of A into e for the
|
||||
// subsequent calculation of the row transformation.
|
||||
e[j] = A[k][j];
|
||||
}
|
||||
if (k < nct) {
|
||||
// Place the transformation in U for subsequent back
|
||||
// multiplication.
|
||||
for (int i = k; i < m; i++) {
|
||||
U[i][k] = A[i][k];
|
||||
}
|
||||
}
|
||||
if (k < nrt) {
|
||||
// Compute the k-th row transformation and place the
|
||||
// k-th super-diagonal in e[k].
|
||||
// Compute 2-norm without under/overflow.
|
||||
e[k] = 0;
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
e[k] = FastMath.hypot(e[k], e[i]);
|
||||
}
|
||||
if (e[k] != 0) {
|
||||
if (e[k + 1] < 0) {
|
||||
e[k] = -e[k];
|
||||
}
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
e[i] /= e[k];
|
||||
}
|
||||
e[k + 1] += 1;
|
||||
}
|
||||
e[k] = -e[k];
|
||||
if (k + 1 < m &&
|
||||
e[k] != 0) {
|
||||
// Apply the transformation.
|
||||
for (int i = k + 1; i < m; i++) {
|
||||
work[i] = 0;
|
||||
}
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
for (int i = k + 1; i < m; i++) {
|
||||
work[i] += e[j] * A[i][j];
|
||||
}
|
||||
}
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
final double t = -e[j] / e[k + 1];
|
||||
for (int i = k + 1; i < m; i++) {
|
||||
A[i][j] += t * work[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Place the transformation in V for subsequent
|
||||
// back multiplication.
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
V[i][k] = e[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
// Set up the final bidiagonal matrix or order p.
|
||||
int p = n;
|
||||
if (nct < n) {
|
||||
singularValues[nct] = A[nct][nct];
|
||||
}
|
||||
if (m < p) {
|
||||
singularValues[p - 1] = 0;
|
||||
}
|
||||
if (nrt + 1 < p) {
|
||||
e[nrt] = A[nrt][p - 1];
|
||||
}
|
||||
e[p - 1] = 0;
|
||||
|
||||
// Generate U.
|
||||
for (int j = nct; j < n; j++) {
|
||||
for (int i = 0; i < m; i++) {
|
||||
U[i][j] = 0;
|
||||
}
|
||||
U[j][j] = 1;
|
||||
}
|
||||
for (int k = nct - 1; k >= 0; k--) {
|
||||
if (singularValues[k] != 0) {
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
double t = 0;
|
||||
for (int i = k; i < m; i++) {
|
||||
t += U[i][k] * U[i][j];
|
||||
}
|
||||
t = -t / U[k][k];
|
||||
for (int i = k; i < m; i++) {
|
||||
U[i][j] += t * U[i][k];
|
||||
}
|
||||
}
|
||||
for (int i = k; i < m; i++) {
|
||||
U[i][k] = -U[i][k];
|
||||
}
|
||||
U[k][k] = 1 + U[k][k];
|
||||
for (int i = 0; i < k - 1; i++) {
|
||||
U[i][k] = 0;
|
||||
}
|
||||
} else {
|
||||
for (int i = 0; i < m; i++) {
|
||||
U[i][k] = 0;
|
||||
}
|
||||
U[k][k] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
// Generate V.
|
||||
for (int k = n - 1; k >= 0; k--) {
|
||||
if (k < nrt &&
|
||||
e[k] != 0) {
|
||||
for (int j = k + 1; j < n; j++) {
|
||||
double t = 0;
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
t += V[i][k] * V[i][j];
|
||||
}
|
||||
t = -t / V[k + 1][k];
|
||||
for (int i = k + 1; i < n; i++) {
|
||||
V[i][j] += t * V[i][k];
|
||||
}
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < n; i++) {
|
||||
V[i][k] = 0;
|
||||
}
|
||||
V[k][k] = 1;
|
||||
}
|
||||
|
||||
// Main iteration loop for the singular values.
|
||||
final int pp = p - 1;
|
||||
int iter = 0;
|
||||
while (p > 0) {
|
||||
int k;
|
||||
int kase;
|
||||
// Here is where a test for too many iterations would go.
|
||||
// This section of the program inspects for
|
||||
// negligible elements in the s and e arrays. On
|
||||
// completion the variables kase and k are set as follows.
|
||||
// kase = 1 if s(p) and e[k-1] are negligible and k<p
|
||||
// kase = 2 if s(k) is negligible and k<p
|
||||
// kase = 3 if e[k-1] is negligible, k<p, and
|
||||
// s(k), ..., s(p) are not negligible (qr step).
|
||||
// kase = 4 if e(p-1) is negligible (convergence).
|
||||
for (k = p - 2; k >= 0; k--) {
|
||||
final double threshold
|
||||
= TINY + EPS * (FastMath.abs(singularValues[k]) +
|
||||
FastMath.abs(singularValues[k + 1]));
|
||||
if (FastMath.abs(e[k]) <= threshold) {
|
||||
e[k] = 0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (k == p - 2) {
|
||||
kase = 4;
|
||||
} else {
|
||||
int ks;
|
||||
for (ks = p - 1; ks >= k; ks--) {
|
||||
if (ks == k) {
|
||||
break;
|
||||
}
|
||||
final double t = (ks != p ? FastMath.abs(e[ks]) : 0) +
|
||||
(ks != k + 1 ? FastMath.abs(e[ks - 1]) : 0);
|
||||
if (FastMath.abs(singularValues[ks]) <= TINY + EPS * t) {
|
||||
singularValues[ks] = 0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if (ks == k) {
|
||||
kase = 3;
|
||||
} else if (ks == p - 1) {
|
||||
kase = 1;
|
||||
} else {
|
||||
kase = 2;
|
||||
k = ks;
|
||||
}
|
||||
}
|
||||
k++;
|
||||
// Perform the task indicated by kase.
|
||||
switch (kase) {
|
||||
// Deflate negligible s(p).
|
||||
case 1: {
|
||||
double f = e[p - 2];
|
||||
e[p - 2] = 0;
|
||||
for (int j = p - 2; j >= k; j--) {
|
||||
double t = FastMath.hypot(singularValues[j], f);
|
||||
final double cs = singularValues[j] / t;
|
||||
final double sn = f / t;
|
||||
singularValues[j] = t;
|
||||
if (j != k) {
|
||||
f = -sn * e[j - 1];
|
||||
e[j - 1] = cs * e[j - 1];
|
||||
}
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
t = cs * V[i][j] + sn * V[i][p - 1];
|
||||
V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1];
|
||||
V[i][j] = t;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
// Split at negligible s(k).
|
||||
case 2: {
|
||||
double f = e[k - 1];
|
||||
e[k - 1] = 0;
|
||||
for (int j = k; j < p; j++) {
|
||||
double t = FastMath.hypot(singularValues[j], f);
|
||||
final double cs = singularValues[j] / t;
|
||||
final double sn = f / t;
|
||||
singularValues[j] = t;
|
||||
f = -sn * e[j];
|
||||
e[j] = cs * e[j];
|
||||
|
||||
for (int i = 0; i < m; i++) {
|
||||
t = cs * U[i][j] + sn * U[i][k - 1];
|
||||
U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1];
|
||||
U[i][j] = t;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
// Perform one qr step.
|
||||
case 3: {
|
||||
// Calculate the shift.
|
||||
final double scale = FastMath.max(FastMath.max(FastMath.max(FastMath.max(
|
||||
FastMath.abs(singularValues[p - 1]), FastMath.abs(singularValues[p - 2])), FastMath.abs(e[p - 2])),
|
||||
FastMath.abs(singularValues[k])), FastMath.abs(e[k]));
|
||||
final double sp = singularValues[p - 1] / scale;
|
||||
final double spm1 = singularValues[p - 2] / scale;
|
||||
final double epm1 = e[p - 2] / scale;
|
||||
final double sk = singularValues[k] / scale;
|
||||
final double ek = e[k] / scale;
|
||||
final double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
|
||||
final double c = (sp * epm1) * (sp * epm1);
|
||||
double shift = 0;
|
||||
if (b != 0 ||
|
||||
c != 0) {
|
||||
shift = FastMath.sqrt(b * b + c);
|
||||
if (b < 0) {
|
||||
shift = -shift;
|
||||
}
|
||||
shift = c / (b + shift);
|
||||
}
|
||||
double f = (sk + sp) * (sk - sp) + shift;
|
||||
double g = sk * ek;
|
||||
// Chase zeros.
|
||||
for (int j = k; j < p - 1; j++) {
|
||||
double t = FastMath.hypot(f, g);
|
||||
double cs = f / t;
|
||||
double sn = g / t;
|
||||
if (j != k) {
|
||||
e[j - 1] = t;
|
||||
}
|
||||
f = cs * singularValues[j] + sn * e[j];
|
||||
e[j] = cs * e[j] - sn * singularValues[j];
|
||||
g = sn * singularValues[j + 1];
|
||||
singularValues[j + 1] = cs * singularValues[j + 1];
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
t = cs * V[i][j] + sn * V[i][j + 1];
|
||||
V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1];
|
||||
V[i][j] = t;
|
||||
}
|
||||
t = FastMath.hypot(f, g);
|
||||
cs = f / t;
|
||||
sn = g / t;
|
||||
singularValues[j] = t;
|
||||
f = cs * e[j] + sn * singularValues[j + 1];
|
||||
singularValues[j + 1] = -sn * e[j] + cs * singularValues[j + 1];
|
||||
g = sn * e[j + 1];
|
||||
e[j + 1] = cs * e[j + 1];
|
||||
if (j < m - 1) {
|
||||
for (int i = 0; i < m; i++) {
|
||||
t = cs * U[i][j] + sn * U[i][j + 1];
|
||||
U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1];
|
||||
U[i][j] = t;
|
||||
}
|
||||
}
|
||||
}
|
||||
e[p - 2] = f;
|
||||
iter = iter + 1;
|
||||
}
|
||||
break;
|
||||
// Convergence.
|
||||
default: {
|
||||
// Make the singular values positive.
|
||||
if (singularValues[k] <= 0) {
|
||||
singularValues[k] = singularValues[k] < 0 ? -singularValues[k] : 0;
|
||||
|
||||
for (int i = 0; i <= pp; i++) {
|
||||
V[i][k] = -V[i][k];
|
||||
}
|
||||
}
|
||||
// Order the singular values.
|
||||
while (k < pp) {
|
||||
if (singularValues[k] >= singularValues[k + 1]) {
|
||||
break;
|
||||
}
|
||||
double t = singularValues[k];
|
||||
singularValues[k] = singularValues[k + 1];
|
||||
singularValues[k + 1] = t;
|
||||
if (k < n - 1) {
|
||||
for (int i = 0; i < n; i++) {
|
||||
t = V[i][k + 1];
|
||||
V[i][k + 1] = V[i][k];
|
||||
V[i][k] = t;
|
||||
}
|
||||
}
|
||||
if (k < m - 1) {
|
||||
for (int i = 0; i < m; i++) {
|
||||
t = U[i][k + 1];
|
||||
U[i][k + 1] = U[i][k];
|
||||
U[i][k] = t;
|
||||
}
|
||||
}
|
||||
k++;
|
||||
}
|
||||
iter = 0;
|
||||
p--;
|
||||
}
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// Set the small value tolerance used to calculate rank and pseudo-inverse
|
||||
tol = FastMath.max(m * singularValues[0] * EPS,
|
||||
FastMath.sqrt(MathUtils.SAFE_MIN));
|
||||
|
||||
if (!transposed) {
|
||||
cachedU = MatrixUtils.createRealMatrix(U);
|
||||
cachedV = MatrixUtils.createRealMatrix(V);
|
||||
} else {
|
||||
cachedU = MatrixUtils.createRealMatrix(V);
|
||||
cachedV = MatrixUtils.createRealMatrix(U);
|
||||
}
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getU() {
|
||||
// return the cached matrix
|
||||
return cachedU;
|
||||
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getUT() {
|
||||
if (cachedUt == null) {
|
||||
cachedUt = getU().transpose();
|
||||
}
|
||||
// return the cached matrix
|
||||
return cachedUt;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getS() {
|
||||
if (cachedS == null) {
|
||||
// cache the matrix for subsequent calls
|
||||
cachedS = MatrixUtils.createRealDiagonalMatrix(singularValues);
|
||||
}
|
||||
return cachedS;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] getSingularValues() {
|
||||
return singularValues.clone();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getV() {
|
||||
// return the cached matrix
|
||||
return cachedV;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getVT() {
|
||||
if (cachedVt == null) {
|
||||
cachedVt = getV().transpose();
|
||||
}
|
||||
// return the cached matrix
|
||||
return cachedVt;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix getCovariance(final double minSingularValue) {
|
||||
// get the number of singular values to consider
|
||||
final int p = singularValues.length;
|
||||
int dimension = 0;
|
||||
while (dimension < p &&
|
||||
singularValues[dimension] >= minSingularValue) {
|
||||
++dimension;
|
||||
}
|
||||
|
||||
if (dimension == 0) {
|
||||
throw new NumberIsTooLargeException(LocalizedFormats.TOO_LARGE_CUTOFF_SINGULAR_VALUE,
|
||||
minSingularValue, singularValues[0], true);
|
||||
}
|
||||
|
||||
final double[][] data = new double[dimension][p];
|
||||
getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() {
|
||||
/** {@inheritDoc} */
|
||||
@Override
|
||||
public void visit(final int row, final int column,
|
||||
final double value) {
|
||||
data[row][column] = value / singularValues[row];
|
||||
}
|
||||
}, 0, dimension - 1, 0, p - 1);
|
||||
|
||||
RealMatrix jv = new Array2DRowRealMatrix(data, false);
|
||||
return jv.transpose().multiply(jv);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double getNorm() {
|
||||
return singularValues[0];
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double getConditionNumber() {
|
||||
return singularValues[0] / singularValues[n - 1];
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes the inverse of the condition number.
|
||||
* In cases of rank deficiency, the {@link #getConditionNumber() condition
|
||||
* number} will become undefined.
|
||||
*
|
||||
* @return the inverse of the condition number.
|
||||
*/
|
||||
public double getInverseConditionNumber() {
|
||||
return singularValues[n - 1] / singularValues[0];
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public int getRank() {
|
||||
int r = 0;
|
||||
for (int i = 0; i < singularValues.length; i++) {
|
||||
if (singularValues[i] > tol) {
|
||||
r++;
|
||||
}
|
||||
}
|
||||
return r;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public DecompositionSolver getSolver() {
|
||||
return new Solver(singularValues, getUT(), getV(), getRank() == m, tol);
|
||||
}
|
||||
|
||||
/** Specialized solver. */
|
||||
private static class Solver implements DecompositionSolver {
|
||||
/** Pseudo-inverse of the initial matrix. */
|
||||
private final RealMatrix pseudoInverse;
|
||||
/** Singularity indicator. */
|
||||
private boolean nonSingular;
|
||||
|
||||
/**
|
||||
* Build a solver from decomposed matrix.
|
||||
*
|
||||
* @param singularValues Singular values.
|
||||
* @param uT U<sup>T</sup> matrix of the decomposition.
|
||||
* @param v V matrix of the decomposition.
|
||||
* @param nonSingular Singularity indicator.
|
||||
* @param tol tolerance for singular values
|
||||
*/
|
||||
private Solver(final double[] singularValues, final RealMatrix uT,
|
||||
final RealMatrix v, final boolean nonSingular, final double tol) {
|
||||
final double[][] suT = uT.getData();
|
||||
for (int i = 0; i < singularValues.length; ++i) {
|
||||
final double a;
|
||||
if (singularValues[i] > tol) {
|
||||
a = 1 / singularValues[i];
|
||||
} else {
|
||||
a = 0;
|
||||
}
|
||||
final double[] suTi = suT[i];
|
||||
for (int j = 0; j < suTi.length; ++j) {
|
||||
suTi[j] *= a;
|
||||
}
|
||||
}
|
||||
pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false));
|
||||
this.nonSingular = nonSingular;
|
||||
}
|
||||
|
||||
/**
|
||||
* Solve the linear equation A × X = B in least square sense.
|
||||
* <p>
|
||||
* The m×n matrix A may not be square, the solution X is such that
|
||||
* ||A × X - B|| is minimal.
|
||||
* </p>
|
||||
* @param b Right-hand side of the equation A × X = B
|
||||
* @return a vector X that minimizes the two norm of A × X - B
|
||||
* @throws org.apache.commons.math.exception.DimensionMismatchException
|
||||
* if the matrices dimensions do not match.
|
||||
*/
|
||||
public RealVector solve(final RealVector b) {
|
||||
return pseudoInverse.operate(b);
|
||||
}
|
||||
|
||||
/**
|
||||
* Solve the linear equation A × X = B in least square sense.
|
||||
* <p>
|
||||
* The m×n matrix A may not be square, the solution X is such that
|
||||
* ||A × X - B|| is minimal.
|
||||
* </p>
|
||||
*
|
||||
* @param b Right-hand side of the equation A × X = B
|
||||
* @return a matrix X that minimizes the two norm of A × X - B
|
||||
* @throws org.apache.commons.math.exception.DimensionMismatchException
|
||||
* if the matrices dimensions do not match.
|
||||
*/
|
||||
public RealMatrix solve(final RealMatrix b) {
|
||||
return pseudoInverse.multiply(b);
|
||||
}
|
||||
|
||||
/**
|
||||
* Check if the decomposed matrix is non-singular.
|
||||
*
|
||||
* @return {@code true} if the decomposed matrix is non-singular.
|
||||
*/
|
||||
public boolean isNonSingular() {
|
||||
return nonSingular;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get the pseudo-inverse of the decomposed matrix.
|
||||
*
|
||||
* @return the inverse matrix.
|
||||
*/
|
||||
public RealMatrix getInverse() {
|
||||
return pseudoInverse;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -27,7 +27,7 @@ import org.junit.Assert;
|
|||
import org.junit.Test;
|
||||
|
||||
|
||||
public class SingularValueDecompositionImplTest {
|
||||
public class SingularValueDecompositionTest {
|
||||
|
||||
private double[][] testSquare = {
|
||||
{ 24.0 / 25.0, 43.0 / 25.0 },
|
||||
|
|
@ -50,7 +50,7 @@ public class SingularValueDecompositionImplTest {
|
|||
final int columns = singularValues.length;
|
||||
Random r = new Random(15338437322523l);
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecompositionImpl(createTestMatrix(r, rows, columns, singularValues));
|
||||
new SingularValueDecomposition(createTestMatrix(r, rows, columns, singularValues));
|
||||
double[] computedSV = svd.getSingularValues();
|
||||
Assert.assertEquals(singularValues.length, computedSV.length);
|
||||
for (int i = 0; i < singularValues.length; ++i) {
|
||||
|
|
@ -65,7 +65,7 @@ public class SingularValueDecompositionImplTest {
|
|||
final int columns = singularValues.length + 2;
|
||||
Random r = new Random(732763225836210l);
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecompositionImpl(createTestMatrix(r, rows, columns, singularValues));
|
||||
new SingularValueDecomposition(createTestMatrix(r, rows, columns, singularValues));
|
||||
double[] computedSV = svd.getSingularValues();
|
||||
Assert.assertEquals(singularValues.length, computedSV.length);
|
||||
for (int i = 0; i < singularValues.length; ++i) {
|
||||
|
|
@ -79,7 +79,7 @@ public class SingularValueDecompositionImplTest {
|
|||
RealMatrix matrix = MatrixUtils.createRealMatrix(testSquare);
|
||||
final int m = matrix.getRowDimension();
|
||||
final int n = matrix.getColumnDimension();
|
||||
SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
|
||||
SingularValueDecomposition svd = new SingularValueDecomposition(matrix);
|
||||
Assert.assertEquals(m, svd.getU().getRowDimension());
|
||||
Assert.assertEquals(m, svd.getU().getColumnDimension());
|
||||
Assert.assertEquals(m, svd.getS().getColumnDimension());
|
||||
|
|
@ -98,7 +98,7 @@ public class SingularValueDecompositionImplTest {
|
|||
{ 9.0 / 2.0, 3.0 / 2.0, 15.0 / 2.0, 5.0 / 2.0 },
|
||||
{ 3.0 / 2.0, 9.0 / 2.0, 5.0 / 2.0, 15.0 / 2.0 }
|
||||
}, false);
|
||||
SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
|
||||
SingularValueDecomposition svd = new SingularValueDecomposition(matrix);
|
||||
Assert.assertEquals(16.0, svd.getSingularValues()[0], 1.0e-14);
|
||||
Assert.assertEquals( 8.0, svd.getSingularValues()[1], 1.0e-14);
|
||||
Assert.assertEquals( 4.0, svd.getSingularValues()[2], 1.0e-14);
|
||||
|
|
@ -135,7 +135,7 @@ public class SingularValueDecompositionImplTest {
|
|||
}
|
||||
|
||||
public void checkAEqualUSVt(final RealMatrix matrix) {
|
||||
SingularValueDecomposition svd = new SingularValueDecompositionImpl(matrix);
|
||||
SingularValueDecomposition svd = new SingularValueDecomposition(matrix);
|
||||
RealMatrix u = svd.getU();
|
||||
RealMatrix s = svd.getS();
|
||||
RealMatrix v = svd.getV();
|
||||
|
|
@ -147,17 +147,17 @@ public class SingularValueDecompositionImplTest {
|
|||
/** test that U is orthogonal */
|
||||
@Test
|
||||
public void testUOrthogonal() {
|
||||
checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getU());
|
||||
checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare)).getU());
|
||||
checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare).transpose()).getU());
|
||||
checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare)).getU());
|
||||
checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare)).getU());
|
||||
checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare).transpose()).getU());
|
||||
}
|
||||
|
||||
/** test that V is orthogonal */
|
||||
@Test
|
||||
public void testVOrthogonal() {
|
||||
checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getV());
|
||||
checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare)).getV());
|
||||
checkOrthogonal(new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare).transpose()).getV());
|
||||
checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare)).getV());
|
||||
checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare)).getV());
|
||||
checkOrthogonal(new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare).transpose()).getV());
|
||||
}
|
||||
|
||||
public void checkOrthogonal(final RealMatrix m) {
|
||||
|
|
@ -171,7 +171,7 @@ public class SingularValueDecompositionImplTest {
|
|||
// together, the actual triplet (U,S,V) is not uniquely defined.
|
||||
public void testMatricesValues1() {
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
|
||||
new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare));
|
||||
RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 3.0 / 5.0, -4.0 / 5.0 },
|
||||
{ 4.0 / 5.0, 3.0 / 5.0 }
|
||||
|
|
@ -224,7 +224,7 @@ public class SingularValueDecompositionImplTest {
|
|||
|
||||
// check values against known references
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testNonSquare));
|
||||
new SingularValueDecomposition(MatrixUtils.createRealMatrix(testNonSquare));
|
||||
RealMatrix u = svd.getU();
|
||||
Assert.assertEquals(0, u.subtract(uRef).getNorm(), normTolerance);
|
||||
RealMatrix s = svd.getS();
|
||||
|
|
@ -244,34 +244,34 @@ public class SingularValueDecompositionImplTest {
|
|||
public void testRank() {
|
||||
double[][] d = { { 1, 1, 1 }, { 0, 0, 0 }, { 1, 2, 3 } };
|
||||
RealMatrix m = new Array2DRowRealMatrix(d);
|
||||
SingularValueDecomposition svd = new SingularValueDecompositionImpl(m);
|
||||
Assert.assertEquals(2, svd.getRank());
|
||||
SingularValueDecomposition svd = new SingularValueDecomposition(m);
|
||||
Assert.assertEquals(2, svd.getRank());
|
||||
}
|
||||
|
||||
|
||||
/** test MATH-583 */
|
||||
@Test
|
||||
public void testStability1() {
|
||||
RealMatrix m = new Array2DRowRealMatrix(201, 201);
|
||||
loadRealMatrix(m,"matrix1.csv");
|
||||
try {
|
||||
new SingularValueDecompositionImpl(m);
|
||||
new SingularValueDecomposition(m);
|
||||
} catch (Exception e) {
|
||||
Assert.fail("Exception whilst constructing SVD");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/** test MATH-327 */
|
||||
@Test
|
||||
public void testStability2() {
|
||||
RealMatrix m = new Array2DRowRealMatrix(7, 168);
|
||||
loadRealMatrix(m,"matrix2.csv");
|
||||
try {
|
||||
new SingularValueDecompositionImpl(m);
|
||||
new SingularValueDecomposition(m);
|
||||
} catch (Throwable e) {
|
||||
Assert.fail("Exception whilst constructing SVD");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
private void loadRealMatrix(RealMatrix m, String resourceName) {
|
||||
try {
|
||||
DataInputStream in = new DataInputStream(getClass().getResourceAsStream(resourceName));
|
||||
|
|
@ -286,25 +286,25 @@ public class SingularValueDecompositionImplTest {
|
|||
row++;
|
||||
}
|
||||
in.close();
|
||||
} catch (IOException e) {}
|
||||
} catch (IOException e) {}
|
||||
}
|
||||
|
||||
|
||||
/** test condition number */
|
||||
@Test
|
||||
public void testConditionNumber() {
|
||||
SingularValueDecompositionImpl svd =
|
||||
new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare));
|
||||
// replace 1.0e-15 with 1.5e-15
|
||||
Assert.assertEquals(3.0, svd.getConditionNumber(), 1.5e-15);
|
||||
}
|
||||
|
||||
@Test
|
||||
public void testInverseConditionNumber() {
|
||||
SingularValueDecompositionImpl svd =
|
||||
new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare));
|
||||
Assert.assertEquals(1.0/3.0, svd.getInverseConditionNumber(), 1.5e-15);
|
||||
}
|
||||
|
||||
|
||||
private RealMatrix createTestMatrix(final Random r, final int rows, final int columns,
|
||||
final double[] singularValues) {
|
||||
final RealMatrix u =
|
||||
|
|
@ -35,7 +35,7 @@ public class SingularValueSolverTest {
|
|||
@Test
|
||||
public void testSolveDimensionErrors() {
|
||||
DecompositionSolver solver =
|
||||
new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getSolver();
|
||||
new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare)).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[3][2]);
|
||||
try {
|
||||
solver.solve(b);
|
||||
|
|
@ -65,7 +65,7 @@ public class SingularValueSolverTest {
|
|||
{ 1.0, 0.0 },
|
||||
{ 0.0, 0.0 }
|
||||
});
|
||||
DecompositionSolver solver = new SingularValueDecompositionImpl(m).getSolver();
|
||||
DecompositionSolver solver = new SingularValueDecomposition(m).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 11, 12 }, { 21, 22 }
|
||||
});
|
||||
|
|
@ -86,7 +86,7 @@ public class SingularValueSolverTest {
|
|||
@Test
|
||||
public void testSolve() {
|
||||
DecompositionSolver solver =
|
||||
new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare)).getSolver();
|
||||
new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare)).getSolver();
|
||||
RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
|
||||
{ 1, 2, 3 }, { 0, -5, 1 }
|
||||
});
|
||||
|
|
@ -119,8 +119,8 @@ public class SingularValueSolverTest {
|
|||
/** test condition number */
|
||||
@Test
|
||||
public void testConditionNumber() {
|
||||
SingularValueDecompositionImpl svd =
|
||||
new SingularValueDecompositionImpl(MatrixUtils.createRealMatrix(testSquare));
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecomposition(MatrixUtils.createRealMatrix(testSquare));
|
||||
// replace 1.0e-15 with 1.5e-15
|
||||
Assert.assertEquals(3.0, svd.getConditionNumber(), 1.5e-15);
|
||||
}
|
||||
|
|
@ -131,7 +131,7 @@ public class SingularValueSolverTest {
|
|||
{ 1.0, 2.0 }, { 1.0, 2.0 }
|
||||
});
|
||||
SingularValueDecomposition svd =
|
||||
new SingularValueDecompositionImpl(rm);
|
||||
new SingularValueDecomposition(rm);
|
||||
RealMatrix recomposed = svd.getU().multiply(svd.getS()).multiply(svd.getVT());
|
||||
Assert.assertEquals(0.0, recomposed.subtract(rm).getNorm(), 2.0e-15);
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue