SingularValueDecomposer.java
/*
* Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
*
* Licensed 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 com.irurueta.algebra;
/**
* Computes Singular Value matrix decomposition, which consists
* on factoring provided input matrix into three factors consisting of 2
* unary matrices and 1 diagonal matrix containing singular values,
* following next expression: A = U * S * V'.
* Where A is provided input matrix of size m-by-n, U is an m-by-n unary
* matrix, S is an n-by-n diagonal matrix containing singular values, and V'
* denotes the transpose/conjugate of V and is an n-by-n unary matrix, for
* m < n.
*/
@SuppressWarnings("DuplicatedCode")
public class SingularValueDecomposer extends Decomposer {
/**
* Constant defining default number of iterations to obtain convergence on
* singular value estimation.
*/
public static final int DEFAULT_MAX_ITERS = 30;
/**
* Constant defining minimum number of iterations allowed to obtain
* convergence on singular value estimation.
*/
public static final int MIN_ITERS = 1;
/**
* Constant defining minimum allowed value as threshold to determine
* whether a singular value is negligible or not.
*/
public static final double MIN_THRESH = 0.0;
/**
* Constant defining machine precision.
*/
public static final double EPS = 1e-12;
/**
* Internal storage of U.
*/
private Matrix u;
/**
* Internal storage of V.
*/
private Matrix v;
/**
* Internal storage of singular values.
*/
private double[] w;
/**
* Contains epsilon value, which indicates an estimation of numerical
* precision given by this machine.
*/
private final double eps;
/**
* Contains threshold used to determine whether a singular value can be
* neglected or not due to numerical precision errors. This can be used to
* determine effective rank of input matrix.
*/
private double tsh;
/**
* Member containing maximum number of iterations to obtain convergence of
* singular values estimation.
* If singular values do not converge on provided maximum number of
* iterations, then a NoConvergenceException will be thrown when calling
* decompose() method.
*/
private int maxIters;
/**
* Constructor of this class.
*/
public SingularValueDecomposer() {
super();
maxIters = DEFAULT_MAX_ITERS;
u = v = null;
w = null;
eps = EPS;
}
/**
* Constructor of this class.
*
* @param maxIters Determines maximum number of iterations to be done when
* decomposing input matrix into singular values so that singular values
* converge properly.
*/
public SingularValueDecomposer(final int maxIters) {
super();
this.maxIters = maxIters;
u = v = null;
w = null;
eps = EPS;
}
/**
* Constructor of this class.
*
* @param inputMatrix Reference to input matrix to be decomposed.
*/
public SingularValueDecomposer(final Matrix inputMatrix) {
super(inputMatrix);
maxIters = DEFAULT_MAX_ITERS;
u = v = null;
w = null;
eps = EPS;
}
/**
* Constructor of this class.
*
* @param inputMatrix Reference to input matrix to be decomposed.
* @param maxIters Determines maximum number of iterations to be done when
* decomposing input matrix into singular value so that singular values
* converge properly.
*/
public SingularValueDecomposer(final Matrix inputMatrix, final int maxIters) {
super(inputMatrix);
this.maxIters = maxIters;
u = v = null;
w = null;
eps = EPS;
}
/**
* Returns decomposer type corresponding to Singular Value decomposition.
*
* @return Decomposer type.
*/
@Override
public DecomposerType getDecomposerType() {
return DecomposerType.SINGULAR_VALUE_DECOMPOSITION;
}
/**
* Sets reference to input matrix to be decomposed.
*
* @param inputMatrix Reference to input matrix to be decomposed.
* @throws LockedException Exception thrown if attempting to call this
* method while this instance remains locked.
*/
@Override
public void setInputMatrix(final Matrix inputMatrix) throws LockedException {
super.setInputMatrix(inputMatrix);
u = v = null;
w = null;
}
/**
* Returns boolean indicating whether decomposition has been computed and
* results can be retrieved.
* Attempting to retrieve decomposition results when not available, will
* probably raise a NotAvailableException
*
* @return Boolean indicating whether decomposition has been computed and
* results can be retrieved.
*/
@Override
public boolean isDecompositionAvailable() {
return u != null || v != null || w != null;
}
/**
* This method computes Singular Value matrix decomposition, which consists
* on factoring provided input matrix into three factors consisting of 2
* unary matrices and 1 diagonal matrix containing singular values,
* following next expression: A = U * S * V'.
* Where A is provided input matrix of size m-by-n, U is an m-by-n unary
* matrix, S is an n-by-n diagonal matrix containing singular values, and V'
* denotes the transpose/conjugate of V and is an n-by-n unary matrix, for
* m < n.
* Note: Factors U, S and V will be accessible once Singular Value
* decomposition has been computed.
* Note: During execution of this method, Singular Value decomposition will
* be available and operations such as retrieving matrix factors, or
* computing rank of matrices among others will be able to be done.
* Attempting to call any of such operations before calling this method will
* raise a NotAvailableException because they require computation of
* SingularValue decomposition first.
*
* @throws NotReadyException Exception thrown if attempting to call this
* method when this instance is not ready (i.e. no input matrix has been
* provided).
* @throws LockedException Exception thrown if this decomposer is already
* locked before calling this method. Notice that this method will actually
* lock this instance while it is being executed.
* @throws DecomposerException Exception thrown if for any reason
* decomposition fails while being executed, like when convergence of
* results cannot be obtained, etc.
*/
@Override
public void decompose() throws NotReadyException, LockedException, DecomposerException {
if (!isReady()) {
throw new NotReadyException();
}
if (isLocked()) {
throw new LockedException();
}
locked = true;
final var m = inputMatrix.getRows();
final var n = inputMatrix.getColumns();
// copy input matrix into U
u = new Matrix(inputMatrix);
w = new double[n];
try {
v = new Matrix(n, n);
} catch (final WrongSizeException ignore) {
// never happens
}
try {
internalDecompose();
reorder();
setNegligibleSingularValueThreshold(0.5 * Math.sqrt(m + n + 1.0) * w[0] * eps);
locked = false;
} catch (final DecomposerException e) {
u = v = null;
w = null;
locked = false;
throw e;
}
}
/**
* Returns maximum number of iterations to be done in order to obtain
* convergence of singular values when computing input matrix Singular
* Value Decomposition.
*
* @return Maximum number of iterations to obtain singular values
* convergence.
*/
public int getMaxIterations() {
return maxIters;
}
/**
* SSets maximum number of iterations to be done in order to obtain
* convergence of singular values when computing input matrix Singular
* Value Decomposition.
* Note: This parameter should rarely be modified because default value
* is usually good enough.
* Note: If convergence of singular values is not achieved within provided
* maximum number of iterations, a DecomposerException will be thrown when
* calling decompose();
*
* @param maxIters Maximum number of iterations to obtain convergence of
* singular values. Provided value must be 1 or greater, otherwise an
* IllegalArgumentException will be thrown.
* @throws LockedException Exception thrown if attempting to call this
* method while this instance remains locked.
* @throws IllegalArgumentException Exception thrown if provided value
* for maxIters is out of valid range of values.
*/
public void setMaxIterations(final int maxIters) throws LockedException {
if (isLocked()) {
throw new LockedException();
}
if (maxIters < MIN_ITERS) {
throw new IllegalArgumentException();
}
this.maxIters = maxIters;
}
/**
* Returns threshold to be used for determining whether a singular value is
* negligible or not.
* This threshold can be used to consider a singular value as zero or not,
* since small singular values might appear in places where they should be
* zero because of rounding errors and machine precision.
* Singular values considered as zero determine aspects such as rank,
* nullability, null-space or range space.
*
* @return Threshold to be used for determining whether a singular value
* is negligible or not.
* @throws NotAvailableException Exception thrown if attempting to call
* this method before computing Singular Value decomposition.
* To avoid this exception call decompose() method first.
*/
public double getNegligibleSingularValueThreshold() throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
return tsh;
}
/**
* Returns a new matrix instance containing the left singular vector (U
* factor) from Singular Value matrix decomposition, which consists
* on decomposing a matrix using the following expression:
* A = U * S * V'.
* Where A is provided input matrix of size m-by-n and U is an m-by-n
* unary matrix for m < n.
*
* @return Matrix instance containing the left singular vectors from a
* Singular Value decomposition.
* @throws NotAvailableException Exception thrown if attempting to call
* this method before computing Singular Value decomposition. To avoid
* this exception call decompose() method first.
* @see #decompose()
*/
public Matrix getU() throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
return u;
}
/**
* Returns a new matrix instance containing the right singular vectors
* (V factor) from Singular Value matrix decomposition, which consists on
* decomposing a matrix using the following expression:
* A = U * S * V',
* Where A is provided input matrix of size m-by-n and V' denotes the
* transpose/conjugate of V, which is an n-by-n unary matrix for m < n.
*
* @return Matrix instance containing the right singular vectors from a
* Singular Value decomposition.
* @throws NotAvailableException Exception thrown if attempting to call
* this method before computing Singular Value decomposition. To avoid
* this exception call decompose() method first.
* @see #decompose()
*/
public Matrix getV() throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
return v;
}
/**
* Returns a new vector instance containing all singular values after
* decomposition.
* Returned vector is equal to the diagonal of S matrix within expression:
* A = U * S * V' where A is provided input matrix and S is a diagonal
* matrix containing singular values on its diagonal.
*
* @return singular values.
* @throws NotAvailableException if decomposition has not yet been computed.
*/
public double[] getSingularValues() throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
return w;
}
/**
* Copies diagonal matrix into provided instance containing all singular
* values on its diagonal after Singular Value matrix decomposition, which
* consists on decomposing a matrix using the following expression:
* A = U * S * V'.
* Where A is provided input matrix of size m-by-n and S is a diagonal
* matrix of size n-by-n for m < n.
*
* @param m matrix instance containing all singular values on its
* diagonal after execution of this method.
* @throws NotAvailableException Exception thrown if attempting to call
* this method before computing Singular Value decomposition. To avoid
* this exception call decompose() method first.
* @throws WrongSizeException if provided matrix does not have size n-by-n.
* @see #decompose()
*/
public void getW(final Matrix m) throws NotAvailableException, WrongSizeException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
if (m.getRows() != w.length || m.getColumns() != w.length) {
throw new WrongSizeException();
}
Matrix.diagonal(w, m);
}
/**
* Returns a new diagonal matrix instance containing all singular values on
* its diagonal after Singular Value matrix decomposition, which consists
* on decomposing a matrix using the following expression: A = U * S * V'.
* Where A is provided input matrix of size m-by-n and S is a diagonal
* matrix of size n-by-n for m < n.
*
* @return Returned matrix instance containing all singular values on its
* diagonal.
* @throws NotAvailableException Exception thrown if attempting to call
* this method before computing Singular Value decomposition. To avoid
* this exception call decompose() method first.
* @see #decompose()
*/
public Matrix getW() throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
// copy S array into a new diagonal matrix
return Matrix.diagonal(w);
}
/**
* Returns the 2-norm of provided input matrix, which is equal to the highest
* singular value found after decomposition. This is also called the Ky Fan
* 1-norm.
* This norm is also equal to the square root of Frobenius norm of the
* squared provided input matrix. In other words:
* sqrt(norm(A' * A, 'fro')) in Matlab notation.
* Where A is provided input matrix and A' is its transpose, and hence
* A' * A can be considered the squared matrix of A, and Frobenius norm
* is defined as the square root of the sum of the squared elements of a
* matrix: sqr(sum(A(:).^2))
*
* @return The 2-norm of provided input matrix, which is equal to the
* highest singular value.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @see #decompose()
*/
public double getNorm2() throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
return w[0];
}
/**
* Returns the condition number of provided input matrix found after
* decomposition.
* The condition number of a matrix measures the sensitivity of the
* solution of a system of linear equations to errors in the data.
* It gives an indication of the accuracy of the results from matrix
* inversion and the solution of a linear system of equations.
* A problem with a low condition number is said to be well-conditioned,
* whereas a problem with a high condition number is said to be
* ill-conditioned.
* The condition number is a property of a matrix and is not related to the
* algorithm or floating point accuracy of a machine to solve a linear
* system of equations or make matrix inversion.
* When solving a linear system of equations (A * X = b), one should think
* of the condition number as being (very roughly) the rate at which the
* solution x will change with respect to a change in b.
* Thus, if the condition number is large, even a small error in b may
* cause a large error in x. On the other hand, if the condition number is
* small then the error in x will not be much bigger than the error in b.
* One way to find the condition number is by using the ration of the
* maximal and minimal singular values of a matrix, which is what this
* method returns.
*
* @return The condition number of provided input matrix.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @see #decompose()
*/
public double getConditionNumber() throws NotAvailableException {
return 1.0 / getReciprocalConditionNumber();
}
/**
* Returns the inverse of the condition number, i.e. 1.0 / condition number.
* Hence, when reciprocal condition number is close to zero, input matrix
* will be ill-conditioned.
* For more information see getConditionNumber()
*
* @return Inverse of the condition number.
* @throws NotAvailableException Exception thrown if attempting to call
* this method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @see #decompose()
* @see #getConditionNumber()
*/
public double getReciprocalConditionNumber() throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var columns = inputMatrix.getColumns();
return (w[0] <= 0.0 || w[columns - 1] <= 0.0) ? 0.0 : w[columns - 1] / w[0];
}
/**
* Returns effective numerical matrix rank.
* By definition rank of a matrix can be found as the number of non-zero
* singular values of such matrix found after decomposition.
* However, rounding error and machine precision may lead to small but non-zero
* singular values in a rank deficient matrix.
* This method tries to cope with such rounding errors by taking into
* account only those non-negligible singular values to determine input
* matrix rank.
* The Rank-nullity theorem states that for a matrix A of size m-by-n then:
* rank(A) + nullity(A) = n
* Where A is input matrix and n is the number of columns of such matrix.
*
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @return Effective numerical matrix rank.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is negative. Returned singular values after decomposition
* are always positive, and hence, provided threshold should be a positive
* value close to zero.
* @see #decompose()
*/
public int getRank(final double singularValueThreshold) throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
if (singularValueThreshold < MIN_THRESH) {
throw new IllegalArgumentException();
}
var r = 0;
for (var aW : w) {
if (aW > singularValueThreshold) {
r++;
}
}
return r;
}
/**
* Returns effective numerical matrix rank.
* By definition rank of a matrix can be found as the number of non-zero
* singular values of such matrix found after decomposition.
* However, rounding error and machine precision may lead to small but non-zero
* singular values in a rank deficient matrix.
* This method tries to cope with such rounding error by taking into account
* only those non-negligible singular values to determine input matrix rank.
* The Rank-nullity theorem states that for a matrix A of size m-by-n then:
* rank(A) + nullity(A) = n
* Where A is input matrix and n is number of columns of such matrix.
* Note: This method makes same actions as int getRank(double)
* except that singular value threshold is automatically computed by taking
* into account input matrix size, maximal singular value and machine
* precision. This threshold is good enough for most situations, and hence
* we discourage setting it manually.
*
* @return Effective numerical matrix rank.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @see #decompose()
*/
public int getRank() throws NotAvailableException {
return getRank(getNegligibleSingularValueThreshold());
}
/**
* Returns effective numerical matrix nullity.
* By definition nullity of a matrix can be found as the number of zero or
* negligible singular values of provided input matrix after decomposition.
* Rounding error and machine precision may lead to small but non-zero
* singular values in a rank deficient matrix.
* This method tries to cope with such rounding error by taking into account
* only those negligible singular values to determine input matrix nullity.
* The Rank-nullity theorem states that for a matrix A of size m-by-n then:
* rank(A) + nullity(A) = n
* Where A is input matrix and n is number of columns of such matrix.
*
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @return Effective numerical matrix nullity.
* @throws NotAvailableException Exception thrown if attempting to call
* this method before computing Singular Value decomposition.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is negative. Returned singular values after decomposition
* are always positive, and hence, provided threshold should be a positive
* value close to zero.
* @see #decompose()
*/
public int getNullity(final double singularValueThreshold) throws NotAvailableException {
final var n = inputMatrix.getColumns();
return n - getRank(singularValueThreshold);
}
/**
* Returns effective numerical matrix nullity.
* By definition nullity of a matrix can be found as the number of zero or
* negligible singular values of provided input matrix after decomposition.
* Rounding error and machine precision may lead to small but non-zero
* singular values in a rank deficient matrix.
* This method tries to cope with such rounding error by taking into account
* only those negligible singular values to determine input matrix nullity.
* The Rank-nullity theorem states that for a matrix A of size m-by-n then:
* rank(A) + nullity(A) = n
* Where A is input matrix and n is number of columns of such matrix.
* Note: This method makes the same actions as int getNullity(double) except
* that singular value threshold is automatically computed by taking into
* account input matrix size, maximal singular value and machine precision.
* This threshold is good enough for most situations, and hence we
* discourage setting it manually.
*
* @return Effective numerical matrix nullity.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @see #decompose()
*/
public int getNullity() throws NotAvailableException {
return getNullity(getNegligibleSingularValueThreshold());
}
/**
* Sets into provided range matrix the Range space of provided input matrix,
* which spans a subspace of dimension equal to the rank of input matrix.
* Range space is equal to the columns of U corresponding to non-negligible
* singular values.
*
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @param range Matrix containing Range space of provided input matrix.
* @throws NotAvailableException Exception thrown if input matrix has rank
* zero or also if attempting to call this method before computing Singular
* Value decomposition. To avoid this exception call decompose() method
* first and make sure that input matrix has non-zero rank.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is negative. Returned singular values after decomposition
* are always positive, and hence, provided threshold should be a positive
* near to zero value.
*/
public void getRange(final double singularValueThreshold, final Matrix range) throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var rank = getRank(singularValueThreshold);
internalGetRange(rank, singularValueThreshold, range);
}
/**
* Sets into provided range matrix the Range space of provided input matrix,
* which spans a subspace of dimension equal to the rank of input matrix.
* Range space is equal to the columns of U corresponding to non-negligible
* singular values.
* This method performs same actions as getRange(double, Matrix) except that
* singular value threshold is automatically computed by taking into account
* input matrix size, maximal singular value and machine precision.
* This threshold is good enough for most situations, and hence we
* discourage setting it manually.
*
* @param range Matrix containing Range space of provided input matrix.
* @throws NotAvailableException Exception thrown if input matrix has rank
* zero or also if attempting to call this method before computing Singular
* Value decomposition. To avoid this exception call decompose() method
* first and make sure that input matrix has non-zero rank.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is negative. Returned singular values after decomposition
* are always positive, and hence, provided threshold should be a positive
* near to zero value.
*/
public void getRange(final Matrix range) throws NotAvailableException {
getRange(getNegligibleSingularValueThreshold(), range);
}
/**
* Returns matrix containing Range space of provided input matrix, which
* spans a subspace of dimension equal to the rank of input matrix.
* Range space is equal to the columns of U corresponding to non-negligible
* singular values.
*
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @return Matrix containing Range space of provided input matrix.
* @throws NotAvailableException Exception thrown if input matrix has rank
* zero or also if attempting to call this method before computing Singular
* Value decomposition. To avoid this exception call decompose() method
* first and make sure that input matrix has non-zero rank.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is negative. Returned singular values after decomposition
* are always positive, and hence, provided threshold should be a positive
* near to zero value.
*/
public Matrix getRange(final double singularValueThreshold) throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var rows = inputMatrix.getRows();
final var rank = getRank(singularValueThreshold);
Matrix out;
try {
out = new Matrix(rows, rank);
} catch (final WrongSizeException e) {
throw new NotAvailableException(e);
}
internalGetRange(rank, singularValueThreshold, out);
return out;
}
/**
* Return matrix containing Range space of provided input matrix, which
* spans a subspace of dimension equal to the rank of input matrix.
* Range space is equal to the columns of U corresponding to non-negligible
* singular values.
* This method performs same actions as Matrix getRange(double) except that
* singular value threshold is automatically computed by taking into account
* input matrix size, maximal singular value and machine precision.
* This threshold is good enough for most situations, and hence we
* discourage setting it manually.
*
* @return Matrix containing Range space of provided input matrix
* @throws NotAvailableException Exception thrown if input matrix has rank
* zero or also if attempting to call this method before computing Singular
* Value decomposition. To avoid this exception call decompose() method
* first and make sure that input matrix has non-zero rank.
*/
public Matrix getRange() throws NotAvailableException {
return getRange(getNegligibleSingularValueThreshold());
}
/**
* Sets into provided matrix null-space of provided input matrix, which spans
* a subspace of dimension equal to the nullity of input matrix. Null-space
* is equal to the columns of V corresponding to negligible singular values.
*
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @param nullspace Matrix containing null-space of provided input matrix.
* @throws NotAvailableException Exception thrown if input matrix has full
* rank, and hence its nullity is zero, or also if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first and make sure that input matrix
* is rank deficient.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is negative. Returned singular values after decomposition
* are always positive, and hence, provided threshold should be a positive
* near to zero value.
*/
public void getNullspace(final double singularValueThreshold, final Matrix nullspace)
throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var nullity = getNullity(singularValueThreshold);
internalGetNullspace(nullity, singularValueThreshold, nullspace);
}
/**
* Sets into provided matrix null-space of provided input matrix, which spans
* a subspace of dimension equal to the nullity of input matrix. Null-space
* is equal to the columns of V corresponding to negligible singular values.
* This method performs same actions as getNullspace(double, Matrix) except
* that singular value threshold is automatically computed by taking into
* account input matrix size, maximal singular value and machine precision.
* This threshold is good enough for most situations, and hence we
* discourage setting it manually.
*
* @param nullspace Matrix containing null-space of provided input matrix.
* @throws NotAvailableException Exception thrown if input matrix has full
* rank, and hence its nullity is zero, or also if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first and make sure that input matrix
* is rank deficient.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is negative. Returned singular values after decomposition
* are always positive, and hence, provided threshold should be a positive
* near to zero value.
*/
public void getNullspace(final Matrix nullspace) throws NotAvailableException {
getNullspace(getNegligibleSingularValueThreshold(), nullspace);
}
/**
* Returns matrix containing null-space of provided input matrix, which spans
* a subspace of dimension equal to the nullity of input matrix. Null-space
* is equal to the columns of V corresponding to negligible singular values.
*
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @return Matrix containing null-space of provided input matrix.
* @throws NotAvailableException Exception thrown if input matrix has full
* rank, and hence its nullity is zero, or also if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first and make sure that input matrix
* is rank deficient.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is negative. Returned singular values after decomposition
* are always positive, and hence, provided threshold should be a positive
* near to zero value.
*/
public Matrix getNullspace(final double singularValueThreshold) throws NotAvailableException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
final var columns = inputMatrix.getColumns();
final var nullity = getNullity(singularValueThreshold);
Matrix out;
try {
out = new Matrix(columns, nullity);
} catch (final WrongSizeException e) {
throw new NotAvailableException(e);
}
internalGetNullspace(nullity, singularValueThreshold, out);
return out;
}
/**
* Returns matrix containing null-space of provided input matrix, which
* spans a subspace of dimension equal to the nullity of input matrix.
* Null-space is equal to the columns of V corresponding to negligible
* singular values.
*
* @return Matrix containing null-space of provided input matrix.
* This method performs same actions as Matrix {@link #getNullspace(double)} except
* that singular value threshold is automatically computed by taking into
* account input matrix size, maximal singular value and machine precision.
* This threshold is good enough for most situations, and hence we
* discourage setting it manually.
* @throws NotAvailableException Exception thrown if input matrix has full
* rank, and hence its nullity is zero, or also if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first and make sure that input matrix
* is rank deficient.
*/
public Matrix getNullspace() throws NotAvailableException {
return getNullspace(getNegligibleSingularValueThreshold());
}
/**
* Solves a linear system of equations of the following form: A * X = B
* using the pseudo-inverse to find the least squares solution.
* Where A is the input matrix provided for Singular Value decomposition,
* X is the solution to the system of equations, and B is the parameters
* matrix.
* Note: This method can be reused for different b matrices without having
* to recompute Singular Value decomposition on the same input matrix.
* Note: Provided b matrix must have the same number of rows as provided
* input matrix A, otherwise a WrongSizeException will be raised.
* Note: In order to execute this method, a Singular Value decomposition
* must be available, otherwise a NotAvailableException will be raised. In
* order to avoid this exception call decompose() method first.
* Note: Provided result matrix will be resized if needed
*
* @param b Parameters matrix that determine a linear system of equations.
* Provided matrix must have the same number of rows as provided input
* matrix for Singular Value decomposition. Besides, each column on
* parameters matrix will represent a new system of equations, whose
* solution will be returned on appropriate column as an output of this
* method.
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @param result Matrix containing least squares solution of linear system
* of equations on each column for each column of provided parameters matrix
* b.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @throws WrongSizeException Exception thrown if provided parameters matrix
* (b) does not have the same number of rows as input matrix being Singular
* Value decomposed.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is lower than minimum allowed value (MIN_THRESH).
* @see #decompose()
*/
public void solve(final Matrix b, final double singularValueThreshold, final Matrix result)
throws NotAvailableException, WrongSizeException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
if (b.getRows() != inputMatrix.getRows()) {
throw new WrongSizeException();
}
if (singularValueThreshold < MIN_THRESH) {
throw new IllegalArgumentException();
}
final var m = inputMatrix.getRows();
final var n = inputMatrix.getColumns();
final var p = b.getColumns();
final var bcol = new double[m];
double[] xx;
// resize result matrix if needed
if (result.getRows() != n || result.getColumns() != p) {
result.resize(n, p);
}
for (var j = 0; j < p; j++) {
for (var i = 0; i < m; i++) {
bcol[i] = b.getElementAt(i, j);
}
xx = solve(bcol, singularValueThreshold);
// set column j of X using values in vector xx
result.setSubmatrix(0, j, n - 1, j, xx);
}
}
/**
* Solves a linear system of equations of the following form: A * X = B
* using the pseudo-inverse to find the least squares solution.
* Where A is the input matrix provided for Singular Value decomposition,
* X is the solution to the system of equations, and B is the parameter
* matrix.
* Note: This method can be reused for different b matrices without having
* to recompute Singular Value decomposition on the same input matrix.
* Note: Provided b matrix must have the same number of rows as provided
* input matrix A, otherwise a WrongSizeException will be raised.
* Note: In order to execute this method, a Singular Value decomposition
* must be available, otherwise a NotAvailableException will be raised. In
* order to avoid this exception call decompose() method first.
* Note: This method performs same actions as Matrix solve(Matrix, double)
* except that singular value threshold is automatically computed by taking
* into account input matrix size, maximal singular value and machine
* precision.
* This threshold is good enough for most situations, and hence we
* discourage setting it manually.
* Note: Provided result matrix will be resized if needed
*
* @param b Parameters matrix that determines a linear system of equations.
* Provided matrix must have the same number of rows as provided input
* matrix for Singular Value decomposition. Besides, each column on
* parameters matrix will represent a new system of equations, whose
* solution will be returned on appropriate column as an output of this method.
* @param result Matrix containing least squares solution of linear system
* of equations on each column for each column of provided parameters matrix
* b.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @throws WrongSizeException Exception thrown if provided parameters matrix
* (b) does not have the same number of rows as input matrix being Singular
* Value decomposed.
* @see #decompose()
*/
public void solve(final Matrix b, final Matrix result) throws NotAvailableException, WrongSizeException {
solve(b, getNegligibleSingularValueThreshold(), result);
}
/**
* Solves a linear system of equations of the following form: A * X = B
* using the pseudo-inverse to find the least squares solution.
* Where A is the input matrix provided for Singular Value decomposition,
* X is the solution to the system of equations, and B is the parameter
* matrix.
* Note: This method can be reused for different b matrices without having
* to recompute Singular Value decomposition on the same input matrix.
* Note: Provided b matrix must have the same number of rows as provided
* input matrix A, otherwise a WrongSizeException will be raised.
* Note: In order to execute this method, a Singular Value decomposition
* must be available, otherwise a NotAvailableException will be raised. In
* order to avoid this exception call decompose() method first.
*
* @param b Parameters matrix that determine a linear system of equations.
* Provided matrix must have the same number of rows as provided input
* matrix for Singular Value decomposition. Besides, each column on
* parameters matrix will represent a new system of equations, whose
* solution will be returned on appropriate column as an output of this
* method.
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @return Matrix containing least squares solution of linear system of
* equations on each column for each column of provided parameters matrix b.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @throws WrongSizeException Exception thrown if provided parameters matrix
* (b) does not have the same number of rows as input matrix being Singular
* Value decomposed.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is lower than minimum allowed value (MIN_THRESH).
* @see #decompose()
*/
public Matrix solve(final Matrix b, final double singularValueThreshold)
throws NotAvailableException, WrongSizeException {
final var n = inputMatrix.getColumns();
final var p = b.getColumns();
final var x = new Matrix(n, p);
solve(b, singularValueThreshold, x);
return x;
}
/**
* Solves a linear system of equations of the following form: A * X = B
* using the pseudo-inverse to find the least squares solution.
* Where A is the input matrix provided for Singular Value decomposition,
* X is the solution to the system of equations, and B is the parameter
* matrix.
* Note: This method can be reused for different b matrices without having
* to recompute Singular Value decomposition on the same input matrix.
* Note: Provided b matrix must have the same number of rows as provided
* input matrix A, otherwise a WrongSizeException will be raised.
* Note: In order to execute this method, a Singular Value decomposition
* must be available, otherwise a NotAvailableException will be raised. In
* order to avoid this exception call decompose() method first.
* Note: This method performs same actions as Matrix solve(Matrix, double)
* except that singular value threshold is automatically computed by taking
* into account input matrix size, maximal singular value and machine
* precision.
* This threshold is good enough for most situations, and hence we
* discourage setting it manually.
*
* @param b Parameters matrix that determines a linear system of equations.
* Provided matrix must have the same number of rows as provided input
* matrix for Singular Value decomposition. Besides, each column on
* parameters matrix will represent a new system of equations, whose
* solution will be returned on appropriate column as an output of this method.
* @return Matrix containing least squares solution of linear system of
* equations on each column for each column of provided parameters matrix b.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing Singular Value decomposition. To avoid this
* exception call decompose() method first.
* @throws WrongSizeException Exception thrown if provided parameters matrix
* (b) does not have the same number of rows as input matrix being Singular
* Value decomposed.
* @see #decompose()
*/
public Matrix solve(final Matrix b) throws NotAvailableException,
WrongSizeException {
return solve(b, getNegligibleSingularValueThreshold());
}
/**
* Solves a linear system of equations of the following form: A * X = B
* using the pseudo-inverse to find the least squares solution.
* Where A i s the input matrix provided for Singular Value decomposition,
* X is the solution to the system of equations, and B is the parameters
* array.
* Note: This method can be reused for different b arrays without having
* to recompute Singular Value decomposition on the same input matrix.
* Note: Provided b array must have the same length as the number of rows
* on provided input matrix A, otherwise a WrongSizeException will be
* raised.
* Note: In order to execute this method, a Singular Value decomposition
* must be available, otherwise a NotAvailableException will be raised. In
* order to avoid this exception call decompose() method first.
*
* @param b Parameters array that determines a linear system of equations.
* Provided array must have the same length as number of rows on provided
* input matrix for Singular Value decomposition.
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @param result Vector where least squares solution of linear system of
* equations for provided parameters array b will be stored.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing SingularValue decomposition.
* To avoid this exception call decompose() method first.
* @throws WrongSizeException Exception thrown if provided parameters array
* (b) does not have the same length as number of rows on input matrix being
* Singular Value decomposed or if provided result array does not have the
* same length as the number of columns on input matrix.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is lower than minimum allowed value (MIN_THRESH).
* @see #decompose()
*/
public void solve(final double[] b, final double singularValueThreshold,
final double[] result) throws NotAvailableException, WrongSizeException {
if (!isDecompositionAvailable()) {
throw new NotAvailableException();
}
if (b.length != inputMatrix.getRows()) {
throw new WrongSizeException();
}
if (singularValueThreshold < MIN_THRESH) {
throw new IllegalArgumentException();
}
final var m = inputMatrix.getRows();
final var n = inputMatrix.getColumns();
if (result.length != n) {
throw new WrongSizeException();
}
double s;
final var tmp = new double[n];
for (var j = 0; j < n; j++) {
s = 0.0;
if (w[j] > singularValueThreshold) {
for (var i = 0; i < m; i++) {
s += u.getElementAt(i, j) * b[i];
}
s /= w[j];
}
tmp[j] = s;
}
for (var j = 0; j < n; j++) {
s = 0.0;
for (var jj = 0; jj < n; jj++) {
s += v.getElementAt(j, jj) * tmp[jj];
}
result[j] = s;
}
}
/**
* Solves a linear system of equations of the following form: A * X = B
* using the pseudo-inverse to find the least squares solution.
* Where A i s the input matrix provided for Singular Value decomposition,
* X is the solution to the system of equations, and B is the parameters
* array.
* Note: This method can be reused for different b arrays without having
* to recompute Singular Value decomposition on the same input matrix.
* Note: Provided b array must have the same length as the number of rows
* on provided input matrix A, otherwise a WrongSizeException will be
* raised.
* Note: In order to execute this method, a Singular Value decomposition
* must be available, otherwise a NotAvailableException will be raised. In
* order to avoid this exception call decompose() method first.
* Note: this method performs same actions as double[] solve(double[],
* double) except that singular value threshold is automatically computed by
* taking into account input matrix size, maximal singular value and machine
* precision.
* This threshold is good enough for most situations, and hence we discourage
* setting it manually.
*
* @param b Parameters array that determines a linear system of equations.
* Provided array must have the same length as number of rows on provided
* input matrix for Singular Value decomposition.
* @param result Vector where least squares solution of linear system of
* equations for provided parameters array b will be stored.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing SingularValue decomposition.
* To avoid this exception call decompose() method first.
* @throws WrongSizeException Exception thrown if provided parameters array
* (b) does not have the same length as number of rows on input matrix being
* Singular Value decomposed or if provided result array does not have the
* same length as the number of columns on input matrix.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is lower than minimum allowed value (MIN_THRESH).
* @see #decompose()
*/
public void solve(final double[] b, final double[] result) throws NotAvailableException, WrongSizeException {
solve(b, getNegligibleSingularValueThreshold(), result);
}
/**
* Solves a linear system of equations of the following form: A * X = B
* using the pseudo-inverse to find the least squares solution.
* Where A i s the input matrix provided for Singular Value decomposition,
* X is the solution to the system of equations, and B is the parameters
* array.
* Note: This method can be reused for different b arrays without having
* to recompute Singular Value decomposition on the same input matrix.
* Note: Provided b array must have the same length as the number of rows
* on provided input matrix A, otherwise a WrongSizeException will be
* raised.
* Note: In order to execute this method, a Singular Value decomposition
* must be available, otherwise a NotAvailableException will be raised. In
* order to avoid this exception call decompose() method first.
*
* @param b Parameters array that determines a linear system of equations.
* Provided array must have the same length as number of rows on provided
* input matrix for Singular Value decomposition.
* @param singularValueThreshold Threshold used to determine whether a
* singular value is negligible or not.
* @return Vector containing least squares solution of linear system of
* equations for provided parameters array b.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing SingularValue decomposition.
* To avoid this exception call decompose() method first.
* @throws WrongSizeException Exception thrown if provided parameters array
* (b) does not have the same length as number of rows on input matrix being
* Singular Value decomposed.
* @throws IllegalArgumentException Exception thrown if provided singular
* value threshold is lower than minimum allowed value (MIN_THRESH).
* @see #decompose()
*/
public double[] solve(final double[] b, final double singularValueThreshold) throws NotAvailableException,
WrongSizeException {
final var n = inputMatrix.getColumns();
final var x = new double[n];
solve(b, singularValueThreshold, x);
return x;
}
/**
* Solves a linear system of equations of the following form: A * X = B
* using the pseudo-inverse to find the least squares solution.
* Where A i s the input matrix provided for Singular Value decomposition,
* X is the solution to the system of equations, and B is the parameters
* array.
* Note: This method can be reused for different b arrays without having
* to recompute Singular Value decomposition on the same input matrix.
* Note: Provided b array must have the same length as the number of rows
* on provided input matrix A, otherwise a WrongSizeException will be
* raised.
* Note: In order to execute this method, a Singular Value decomposition
* must be available, otherwise a NotAvailableException will be raised. In
* order to avoid this exception call decompose() method first.
* Note: this method performs same actions as double[] solve(double[],
* double) except that singular value threshold is automatically computed by
* taking into account input matrix size, maximal singular value and machine
* precision.
* This threshold is good enough for most situations, and hence we discourage
* setting it manually.
*
* @param b Parameters array that determines a linear system of equations.
* Provided array must have the same length as number of rows on provided
* input matrix for Singular Value decomposition.
* @return Vector containing least squares solution of linear system of
* equations for provided parameters array b.
* @throws NotAvailableException Exception thrown if attempting to call this
* method before computing SingularValue decomposition.
* To avoid this exception call decompose() method first.
* @throws WrongSizeException Exception thrown if provided parameters array
* (b) does not have the same length as number of rows on input matrix being
* Singular Value decomposed.
* @see #decompose()
*/
public double[] solve(final double[] b) throws NotAvailableException, WrongSizeException {
return solve(b, getNegligibleSingularValueThreshold());
}
/**
* This method is called internally by decompose(), and actually computes
* Singular Value Decomposition.
* However, algorithm implemented in this algorithm does not ensure that
* singular values are ordered from maximal to minimal, and hence reorder()
* method is called next within decompose() as well.
*
* @throws NoConvergenceException Exception thrown if singular value
* estimation does not converge within provided number of maximum
* iterations.
*/
@SuppressWarnings("DuplicatedCode")
private void internalDecompose() throws NoConvergenceException {
final var m = inputMatrix.getRows();
final var n = inputMatrix.getColumns();
boolean flag;
int i;
int its;
int j;
int jj;
int k;
var l = 0;
var nm = 0;
double anorm;
double c;
double f;
double g;
double h;
double s;
double scale;
double x;
double y;
double z;
var rv1 = new double[n];
// Householder reduction to bi-diagonal form
g = scale = anorm = 0.0;
for (i = 0; i < n; i++) {
l = i + 2;
rv1[i] = scale * g;
g = s = scale = 0.0;
if (i < m) {
for (k = i; k < m; k++) {
scale += Math.abs(u.getElementAt(k, i));
}
if (scale != 0.0) {
for (k = i; k < m; k++) {
u.setElementAt(k, i, u.getElementAt(k, i) / scale);
s += Math.pow(u.getElementAt(k, i), 2.0);
}
f = u.getElementAt(i, i);
g = -sign(Math.sqrt(s), f);
h = f * g - s;
u.setElementAt(i, i, f - g);
for (j = l - 1; j < n; j++) {
for (s = 0.0, k = i; k < m; k++) {
s += u.getElementAt(k, i) * u.getElementAt(k, j);
}
f = s / h;
for (k = i; k < m; k++) {
u.setElementAt(k, j, u.getElementAt(k, j) + f * u.getElementAt(k, i));
}
}
for (k = i; k < m; k++) {
u.setElementAt(k, i, u.getElementAt(k, i) * scale);
}
}
}
w[i] = scale * g;
g = s = scale = 0.0;
if (i + 1 <= m && i + 1 != n) {
for (k = l - 1; k < n; k++) {
scale += Math.abs(u.getElementAt(i, k));
}
if (scale != 0.0) {
for (k = l - 1; k < n; k++) {
u.setElementAt(i, k, u.getElementAt(i, k) / scale);
s += Math.pow(u.getElementAt(i, k), 2.0);
}
f = u.getElementAt(i, l - 1);
g = -sign(Math.sqrt(s), f);
h = f * g - s;
u.setElementAt(i, l - 1, f - g);
for (k = l - 1; k < n; k++) {
rv1[k] = u.getElementAt(i, k) / h;
}
for (j = l - 1; j < m; j++) {
for (s = 0.0, k = l - 1; k < n; k++) {
s += u.getElementAt(j, k) * u.getElementAt(i, k);
}
for (k = l - 1; k < n; k++) {
u.setElementAt(j, k, u.getElementAt(j, k) + s * rv1[k]);
}
}
for (k = l - 1; k < n; k++) {
u.setElementAt(i, k, u.getElementAt(i, k) * scale);
}
}
}
anorm = Math.max(anorm, Math.abs(w[i]) + Math.abs(rv1[i]));
}
// Accumulation of right-hand transformations
for (i = n - 1; i >= 0; i--) {
if (i < (n - 1)) {
if (g != 0.0) {
// Double division to avoid possible underflow.
for (j = l; j < n; j++) {
v.setElementAt(j, i, u.getElementAt(i, j) / u.getElementAt(i, l) / g);
}
for (j = l; j < n; j++) {
for (s = 0.0, k = l; k < n; k++) {
s += u.getElementAt(i, k) * v.getElementAt(k, j);
}
for (k = l; k < n; k++) {
v.setElementAt(k, j, v.getElementAt(k, j) + s * v.getElementAt(k, i));
}
}
}
for (j = l; j < n; j++) {
v.setElementAt(i, j, 0.0);
v.setElementAt(j, i, 0.0);
}
}
v.setElementAt(i, i, 1.0);
g = rv1[i];
l = i;
}
// Accumulation of left-hand transformations
for (i = Math.min(m, n) - 1; i >= 0; i--) {
l = i + 1;
g = w[i];
for (j = l; j < n; j++) {
u.setElementAt(i, j, 0.0);
}
if (g != 0.0) {
g = 1.0 / g;
for (j = l; j < n; j++) {
for (s = 0.0, k = l; k < m; k++) {
s += u.getElementAt(k, i) * u.getElementAt(k, j);
}
f = (s / u.getElementAt(i, i)) * g;
for (k = i; k < m; k++) {
u.setElementAt(k, j, u.getElementAt(k, j) + f * u.getElementAt(k, i));
}
}
for (j = i; j < m; j++) {
u.setElementAt(j, i, u.getElementAt(j, i) * g);
}
} else {
for (j = i; j < m; j++) {
u.setElementAt(j, i, 0.0);
}
}
u.setElementAt(i, i, u.getElementAt(i, i) + 1.0);
}
// Diagonalization of the bi-diagonal form: Loop over singular values and
// over allowed iterations.
for (k = n - 1; k >= 0; k--) {
for (its = 0; its < maxIters; its++) {
flag = true;
// Test for splitting
// Note that rrv1[0] is always zero
for (l = k; l >= 0; l--) {
nm = l - 1;
if (l == 0 || Math.abs(rv1[l]) <= eps * anorm) {
flag = false;
}
if (!flag || Math.abs(w[nm]) <= eps * anorm) {
break;
}
}
// Cancellation of rv1[0] if l > 1
if (flag) {
c = 0.0;
s = 1.0;
for (i = l; i < k + 1; i++) {
f = s * rv1[i];
rv1[i] = c * rv1[i];
if (Math.abs(f) <= eps * anorm) {
break;
}
g = w[i];
h = pythag(f, g);
w[i] = h;
if (h != 0.0) {
h = 1.0 / h;
} else {
h = Double.MAX_VALUE;
}
c = g * h;
s = -f * h;
for (j = 0; j < m; j++) {
y = u.getElementAt(j, nm);
z = u.getElementAt(j, i);
u.setElementAt(j, nm, y * c + z * s);
u.setElementAt(j, i, z * c - y * s);
}
}
}
z = w[k];
// Convergence.
if (l == k) {
// Singular value is made non-negative
if (z < 0.0) {
w[k] = -z;
for (j = 0; j < n; j++) {
v.setElementAt(j, k, -v.getElementAt(j, k));
}
}
break;
}
if (its == maxIters - 1) {
throw new NoConvergenceException();
}
// Shift from bottom 2-by-2 minor.
x = w[l];
nm = k - 1;
y = w[nm];
g = rv1[nm];
h = rv1[k];
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
g = pythag(f, 1.0);
f = ((x - z) * (x + z) + h * ((y / (f + sign(g, f))) - h)) / x;
c = s = 1.0;
// Next QR transformation
for (j = l; j <= nm; j++) {
i = j + 1;
g = rv1[i];
y = w[i];
h = s * g;
g = c * g;
z = pythag(f, h);
rv1[j] = z;
if (z != 0.0) {
c = f / z;
s = h / z;
} else {
c = Math.signum(f) * Double.MAX_VALUE;
s = Math.signum(h) * Double.MAX_VALUE;
}
f = x * c + g * s;
g = g * c - x * s;
h = y * s;
y *= c;
for (jj = 0; jj < n; jj++) {
x = v.getElementAt(jj, j);
z = v.getElementAt(jj, i);
v.setElementAt(jj, j, x * c + z * s);
v.setElementAt(jj, i, z * c - x * s);
}
z = pythag(f, h);
// Rotation can be arbitrary if z = 0
w[j] = z;
if (z != 0.0) {
z = 1.0 / z;
c = f * z;
s = h * z;
}
f = c * g + s * y;
x = c * y - s * g;
for (jj = 0; jj < m; jj++) {
y = u.getElementAt(jj, j);
z = u.getElementAt(jj, i);
u.setElementAt(jj, j, y * c + z * s);
u.setElementAt(jj, i, z * c - y * s);
}
}
rv1[l] = 0.0;
rv1[k] = f;
w[k] = x;
}
}
}
/**
* Reorders singular values from maximal to minimal, and also reorders
* columns and rows of U and V to ensure that Singular Value Decomposition
* still remains valid.
*/
private void reorder() {
final var m = inputMatrix.getRows();
final var n = inputMatrix.getColumns();
int i;
int j;
int k;
int s;
var inc = 1;
double sw;
final var su = new double[m];
final var sv = new double[n];
do {
inc *= 3;
inc++;
} while (inc <= n);
do {
inc /= 3;
for (i = inc; i < n; i++) {
sw = w[i];
for (k = 0; k < m; k++) {
su[k] = u.getElementAt(k, i);
}
for (k = 0; k < n; k++) {
sv[k] = v.getElementAt(k, i);
}
j = i;
while (w[j - inc] < sw) {
w[j] = w[j - inc];
for (k = 0; k < m; k++) {
u.setElementAt(k, j, u.getElementAt(k, j - inc));
}
for (k = 0; k < n; k++) {
v.setElementAt(k, j, v.getElementAt(k, j - inc));
}
j -= inc;
if (j < inc) {
break;
}
}
w[j] = sw;
for (k = 0; k < m; k++) {
u.setElementAt(k, j, su[k]);
}
for (k = 0; k < n; k++) {
v.setElementAt(k, j, sv[k]);
}
}
} while (inc > 1);
for (k = 0; k < n; k++) {
s = 0;
for (i = 0; i < m; i++) {
if (u.getElementAt(i, k) < 0.0) {
s++;
}
}
for (j = 0; j < n; j++) {
if (v.getElementAt(j, k) < 0.0) {
s++;
}
}
if (s > (m + n) / 2) {
for (i = 0; i < m; i++) {
u.setElementAt(i, k, -u.getElementAt(i, k));
}
for (j = 0; j < n; j++) {
v.setElementAt(j, k, -v.getElementAt(j, k));
}
}
}
}
/**
* Sets threshold to be used to determine whether a singular value is
* negligible or not.
* This threshold can be used to consider a singular value as zero or not,
* since small singular values might appear in places where they should be
* zero because of rounding errors and machine precision.
* Singular values considered as zero determine aspects such as rank,
* nullability, null-space or range space.
*
* @param threshold Threshold to be used to determine whether a singular
* value is negligible or not.
*/
private void setNegligibleSingularValueThreshold(final double threshold) {
tsh = threshold;
}
/**
* Computes norm of a vector of 2 components 'a' and 'b' as
* sqrt(pow(a, 2.0) + pow(b, 2.0)) without destructive underflow or
* overflow, that is when a or b are close to maximum or minimum values
* allowed by machine precision, computing the previous expression might
* lead to highly inaccurate results.
* This method implements previous expression to avoid this effect as
* much as possible and increase accuracy.
*
* @param a 1st value
* @param b 2nd value
* @return Norm of (a, b).
*/
private double pythag(final double a, final double b) {
final var absa = Math.abs(a);
final var absb = Math.abs(b);
if (absa > absb) {
return absa * Math.sqrt(1.0 + (absb / absa) * (absb / absa));
} else {
return (absb == 0.0 ? 0.0 : absb * Math.sqrt(1.0 + (absa / absb) * (absa / absb)));
}
}
/**
* Returns a or -a depending on b sign. If b is positive, this method
* returns "a", otherwise it returns -a
*
* @param a 1st value
* @param b 2nd value
* @return a or -a depending on b sign.
*/
private double sign(final double a, final double b) {
if (b >= 0.0) {
return a >= 0.0 ? a : -a;
} else {
return a >= 0.0 ? -a : a;
}
}
/**
* Internal method to copy range space vector values into provided matrix.
* Provided matrix will be resized if needed
*
* @param rank Rank of range space
* @param singularValueThreshold Threshold to determine whether a singular
* value is null
* @param range Matrix where range space vector values are stored.
*/
private void internalGetRange(final int rank, final double singularValueThreshold, final Matrix range) {
final var rows = inputMatrix.getRows();
final var columns = inputMatrix.getColumns();
if (range.getRows() != rows || range.getColumns() != rank) {
try {
range.resize(rows, rank);
} catch (final WrongSizeException ignore) {
// never happens
}
}
var nr = 0;
for (var j = 0; j < columns; j++) {
if (w[j] > singularValueThreshold) {
// copy column j of U matrix into column nr of out matrix
range.setSubmatrix(0, nr, rows - 1, nr, u, 0, j,
rows - 1, j);
nr++;
}
}
}
/**
* Internal method to copy null-space vector values into provided matrix.
* Provided matrix will be resized if needed
*
* @param nullity Nullity of null-space
* @param singularValueThreshold Threshold to determine whether a singular
* value is null
* @param nullspace Matrix where null-space vector values are stored.
*/
private void internalGetNullspace(final int nullity, final double singularValueThreshold,
final Matrix nullspace) {
final int columns = inputMatrix.getColumns();
if (nullspace.getRows() != columns || nullspace.getColumns() != nullity) {
try {
nullspace.resize(columns, nullity);
} catch (final WrongSizeException ignore) {
// never happens
}
}
var nn = 0;
for (var j = 0; j < columns; j++) {
if (w[j] <= singularValueThreshold) {
// copy column j of U matrix into column nn of out matrix
nullspace.setSubmatrix(0, nn, columns - 1, nn, v, 0, j,
columns - 1, j);
nn++;
}
}
}
}