Java Code Examples for org.apache.commons.math3.linear.DecompositionSolver#isNonSingular()

/**
 * @param data dense matrix represented in row-major form
 * @return solver for the system Ax = b
 */
static Solver getSolver(double[][] data) {
  if (data == null) {
    return null;
  }
  RealMatrix M = new Array2DRowRealMatrix(data, false);
  double infNorm = M.getNorm();
  double singularityThreshold = infNorm * SINGULARITY_THRESHOLD_RATIO;
  RRQRDecomposition decomposition = new RRQRDecomposition(M, singularityThreshold);
  DecompositionSolver solver = decomposition.getSolver();
  if (solver.isNonSingular()) {
    return new Solver(solver);
  }
  // Otherwise try to report apparent rank
  int apparentRank = decomposition.getRank(0.01); // Better value?
  log.warn("{} x {} matrix is near-singular (threshold {}). Add more data or decrease the " +
           "number of features, to <= about {}",
           M.getRowDimension(), 
           M.getColumnDimension(),
           singularityThreshold,
           apparentRank);
  throw new SingularMatrixSolverException(apparentRank, "Apparent rank: " + apparentRank);
}
 

@Override
public Solver getSolver(RealMatrix M) {
  if (M == null) {
    return null;
  }
  RRQRDecomposition decomposition = new RRQRDecomposition(M, SINGULARITY_THRESHOLD);
  DecompositionSolver solver = decomposition.getSolver();
  if (solver.isNonSingular()) {
    return new CommonsMathSolver(solver);
  }
  // Otherwise try to report apparent rank
  int apparentRank = decomposition.getRank(0.01); // Better value?
  log.warn("{} x {} matrix is near-singular (threshold {}). Add more data or decrease the value of model.features, " +
           "to <= about {}",
           M.getRowDimension(), 
           M.getColumnDimension(), 
           SINGULARITY_THRESHOLD,
           apparentRank);
  throw new SingularMatrixSolverException(apparentRank, "Apparent rank: " + apparentRank);
}
 

/**
 * pdf(x, x_hat) = exp(-0.5 * (x-x_hat) * inv(Σ) * (x-x_hat)T) / ( 2π^0.5d * det(Σ)^0.5)
 * 
 * @return value of probabilistic density function
 * @link https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function
 */
public static double pdf(@Nonnull final RealVector x, @Nonnull final RealVector x_hat,
        @Nonnull final RealMatrix sigma) {
    final int dim = x.getDimension();
    Preconditions.checkArgument(x_hat.getDimension() == dim,
        "|x| != |x_hat|, |x|=" + dim + ", |x_hat|=" + x_hat.getDimension());
    Preconditions.checkArgument(sigma.getRowDimension() == dim,
        "|x| != |sigma|, |x|=" + dim + ", |sigma|=" + sigma.getRowDimension());
    Preconditions.checkArgument(sigma.isSquare(), "Sigma is not square matrix");

    LUDecomposition LU = new LUDecomposition(sigma);
    final double detSigma = LU.getDeterminant();
    double denominator = Math.pow(2.d * Math.PI, 0.5d * dim) * Math.pow(detSigma, 0.5d);
    if (denominator == 0.d) { // avoid divide by zero
        return 0.d;
    }

    final RealMatrix invSigma;
    DecompositionSolver solver = LU.getSolver();
    if (solver.isNonSingular() == false) {
        SingularValueDecomposition svd = new SingularValueDecomposition(sigma);
        invSigma = svd.getSolver().getInverse(); // least square solution
    } else {
        invSigma = solver.getInverse();
    }
    //EigenDecomposition eigen = new EigenDecomposition(sigma);
    //double detSigma = eigen.getDeterminant();
    //RealMatrix invSigma = eigen.getSolver().getInverse();

    RealVector diff = x.subtract(x_hat);
    RealVector premultiplied = invSigma.preMultiply(diff);
    double sum = premultiplied.dotProduct(diff);
    double numerator = Math.exp(-0.5d * sum);

    return numerator / denominator;
}
 

@Nonnull
public static RealMatrix inverse(@Nonnull final RealMatrix m, final boolean exact)
        throws SingularMatrixException {
    LUDecomposition LU = new LUDecomposition(m);
    DecompositionSolver solver = LU.getSolver();
    final RealMatrix inv;
    if (exact || solver.isNonSingular()) {
        inv = solver.getInverse();
    } else {
        SingularValueDecomposition SVD = new SingularValueDecomposition(m);
        inv = SVD.getSolver().getInverse();
    }
    return inv;
}
 

/**
 * L = A x R
 *
 * @return a matrix A that minimizes A x R - L
 */
@Nonnull
public static RealMatrix solve(@Nonnull final RealMatrix L, @Nonnull final RealMatrix R,
        final boolean exact) throws SingularMatrixException {
    LUDecomposition LU = new LUDecomposition(R);
    DecompositionSolver solver = LU.getSolver();
    final RealMatrix A;
    if (exact || solver.isNonSingular()) {
        A = LU.getSolver().solve(L);
    } else {
        SingularValueDecomposition SVD = new SingularValueDecomposition(R);
        A = SVD.getSolver().solve(L);
    }
    return A;
}
 

/**
 * Calculates pseudo inverse of a matrix using QR decomposition
 * @param arr the array to invert
 * @return the pseudo inverted matrix
 */
public static INDArray pinvert(INDArray arr, boolean inPlace) {

    // TODO : do it natively instead of relying on commons-maths

    RealMatrix realMatrix = CheckUtil.convertToApacheMatrix(arr);
    QRDecomposition decomposition = new QRDecomposition(realMatrix, 0);
    DecompositionSolver solver = decomposition.getSolver();

    if (!solver.isNonSingular()) {
        throw new IllegalArgumentException("invalid array: must be singular matrix");
    }

    RealMatrix pinvRM = solver.getInverse();

    INDArray pseudoInverse = CheckUtil.convertFromApacheMatrix(pinvRM);

    if (inPlace)
        arr.assign(pseudoInverse);
    return pseudoInverse;

}
 

/**
 * Calculates pseudo inverse of a matrix using QR decomposition
 * @param arr the array to invert
 * @return the pseudo inverted matrix
 */
public static INDArray pinvert(INDArray arr, boolean inPlace) {

    // TODO : do it natively instead of relying on commons-maths

    RealMatrix realMatrix = CheckUtil.convertToApacheMatrix(arr);
    QRDecomposition decomposition = new QRDecomposition(realMatrix, 0);
    DecompositionSolver solver = decomposition.getSolver();

    if (!solver.isNonSingular()) {
        throw new IllegalArgumentException("invalid array: must be singular matrix");
    }

    RealMatrix pinvRM = solver.getInverse();

    INDArray pseudoInverse = CheckUtil.convertFromApacheMatrix(pinvRM, arr.dataType());

    if (inPlace)
        arr.assign(pseudoInverse);
    return pseudoInverse;

}
 

@Override
public boolean isNonSingular(RealMatrix M) {
  QRDecomposition decomposition = new RRQRDecomposition(M, SINGULARITY_THRESHOLD);
  DecompositionSolver solver = decomposition.getSolver();
  return solver.isNonSingular();
}