org.apache.commons.math3.linear.LUDecomposition Java Examples

/**
 * Inverts a matrix
 * @param arr the array to invert
 * @param inPlace Whether to store the result in {@code arr}
 * @return the inverted matrix
 */
public static INDArray invert(INDArray arr, boolean inPlace) {
    if (!arr.isSquare()) {
        throw new IllegalArgumentException("invalid array: must be square matrix");
    }

    //FIX ME: Please
   /* int[] IPIV = new int[arr.length() + 1];
    int LWORK = arr.length() * arr.length();
    INDArray WORK = Nd4j.create(new double[LWORK]);
    INDArray inverse = inPlace ? arr : arr.dup();
    Nd4j.getBlasWrapper().lapack().getrf(arr);
    Nd4j.getBlasWrapper().lapack().getri(arr.size(0),inverse,arr.size(0),IPIV,WORK,LWORK,0);*/

    RealMatrix rm = CheckUtil.convertToApacheMatrix(arr);
    RealMatrix rmInverse = new LUDecomposition(rm).getSolver().getInverse();


    INDArray inverse = CheckUtil.convertFromApacheMatrix(rmInverse);
    if (inPlace)
        arr.assign(inverse);
    return inverse;

}
 

/**
 * Function to perform LU decomposition on a given matrix.
 * 
 * @param in matrix object
 * @return array of matrix blocks
 */
private static MatrixBlock[] computeLU(MatrixBlock in) {
	if ( in.getNumRows() != in.getNumColumns() ) {
		throw new DMLRuntimeException("LU Decomposition can only be done on a square matrix. Input matrix is rectangular (rows=" + in.getNumRows() + ", cols="+ in.getNumColumns() +")");
	}
	
	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);
	
	// Perform LUP decomposition
	LUDecomposition ludecompose = new LUDecomposition(matrixInput);
	RealMatrix P = ludecompose.getP();
	RealMatrix L = ludecompose.getL();
	RealMatrix U = ludecompose.getU();
	
	// Read the results into native format
	MatrixBlock mbP = DataConverter.convertToMatrixBlock(P.getData());
	MatrixBlock mbL = DataConverter.convertToMatrixBlock(L.getData());
	MatrixBlock mbU = DataConverter.convertToMatrixBlock(U.getData());

	return new MatrixBlock[] { mbP, mbL, mbU };
}
 

/**
 * Calculates the LUP-decomposition of a square matrix. The
 * LUP-decomposition of a matrix A consists of three matrices L, U and P
 * that satisfy: P×A = L×U. L is lower triangular (with unit diagonal
 * terms), U is upper triangular and P is a permutation matrix. All matrices
 * are m×m.
 *
 * @param a Given matrix.
 * @return Result P/L/U arrays.
 */
public static Array[] lu(Array a) {
    Array Pa = Array.factory(DataType.DOUBLE, a.getShape());
    Array La = Array.factory(DataType.DOUBLE, a.getShape());
    Array Ua = Array.factory(DataType.DOUBLE, a.getShape());
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray_Double(a);
    RealMatrix matrix = new Array2DRowRealMatrix(aa, false);
    LUDecomposition decomposition = new LUDecomposition(matrix);
    RealMatrix P = decomposition.getP();
    RealMatrix L = decomposition.getL();
    RealMatrix U = decomposition.getU();
    int n = L.getColumnDimension();
    int m = L.getRowDimension();
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            Pa.setDouble(i * n + j, P.getEntry(i, j));
            La.setDouble(i * n + j, L.getEntry(i, j));
            Ua.setDouble(i * n + j, U.getEntry(i, j));
        }
    }

    return new Array[]{Pa, La, Ua};
}
 

/**
 * @param mu1 mean vector of the first normal distribution
 * @param sigma1 covariance matrix of the first normal distribution
 * @param mu2 mean vector of the second normal distribution
 * @param sigma2 covariance matrix of the second normal distribution
 * @return the Hellinger distance between two multivariate normal distributions
 * @link https://en.wikipedia.org/wiki/Hellinger_distance#Examples
 */
public static double hellingerDistance(@Nonnull final RealVector mu1,
        @Nonnull final RealMatrix sigma1, @Nonnull final RealVector mu2,
        @Nonnull final RealMatrix sigma2) {
    RealVector muSub = mu1.subtract(mu2);
    RealMatrix sigmaMean = sigma1.add(sigma2).scalarMultiply(0.5d);
    LUDecomposition LUsigmaMean = new LUDecomposition(sigmaMean);
    double denominator = Math.sqrt(LUsigmaMean.getDeterminant());
    if (denominator == 0.d) {
        return 1.d; // avoid divide by zero
    }
    RealMatrix sigmaMeanInv = LUsigmaMean.getSolver().getInverse(); // has inverse iff det != 0
    double sigma1Det = MatrixUtils.det(sigma1);
    double sigma2Det = MatrixUtils.det(sigma2);

    double numerator = Math.pow(sigma1Det, 0.25d) * Math.pow(sigma2Det, 0.25d)
            * Math.exp(-0.125d * sigmaMeanInv.preMultiply(muSub).dotProduct(muSub));
    return 1.d - numerator / denominator;
}
 

public double[] predict(LearningInstance instance) {
    RealMatrix A = variableSpace.getMatrixVarByName(LinearUCBKey.A.get());
    RealVector B = variableSpace.getScalarVarByName(LinearUCBKey.B.get());
    RealMatrix invA = new LUDecomposition(A).getSolver().getInverse();
    RealVector theta = invA.operate(B);
    RealVector x = extractDenseVector(theta.getDimension(), instance);
    double bound = Math.sqrt(x.dotProduct(invA.operate(x)));
    double mean = x.dotProduct(theta);
    double pred = mean + alpha * bound;
    if (Double.isNaN(pred)) {
        logger.error("Prediction is NaN, model parameter A probably goes wrong.");
        pred = 0.0;
    }
    double[] preds = new double[1];
    preds[0] = pred;
    return preds;
}
 

/**
 * Returns the log posterior PDF of a particular point xi, to belong to this
 * cluster.
 * 
 * @param xi    The point for which we want to estimate the PDF.
 * @return      The log posterior PDF
 */
@Override
public double posteriorLogPdf(Point xi) {
    RealVector x_mu = xi.data.subtract(mean);

    if(cache_covariance_determinant==null || cache_covariance_inverse==null) {
        LUDecomposition lud = new LUDecomposition(covariance);
        cache_covariance_determinant = lud.getDeterminant();
        cache_covariance_inverse = lud.getSolver().getInverse();
        lud =null;
    }
    Double determinant=cache_covariance_determinant;
    RealMatrix invCovariance=cache_covariance_inverse;

    double x_muInvSx_muT = (invCovariance.preMultiply(x_mu)).dotProduct(x_mu);

    double normConst = 1.0/( Math.pow(2*Math.PI, dimensionality/2.0) * Math.pow(determinant, 0.5) );


    //double pdf = Math.exp(-0.5 * x_muInvSx_muT)*normConst;
    double logPdf = -0.5 * x_muInvSx_muT + Math.log(normConst);
    return logPdf;
}
 

public MultivariateTDistribution(RealVector mean, RealMatrix covarianceMatrix, double degreesOfFreedom) {
    this.mean = mean;
    if (mean.getDimension() > 1) {
        this.precisionMatrix = MatrixUtils.blockInverse(covarianceMatrix, (-1 + covarianceMatrix.getColumnDimension()) / 2);
    } else {
        this.precisionMatrix = MatrixUtils.createRealIdentityMatrix(1).scalarMultiply(1. / covarianceMatrix.getEntry(0, 0));
    }
    this.dof = degreesOfFreedom;

    this.D = mean.getDimension();

    double determinant = new LUDecomposition(covarianceMatrix).getDeterminant();

    this.multiplier = Math.exp(Gamma.logGamma(0.5 * (D + dof)) - Gamma.logGamma(0.5 * dof)) /
            Math.pow(Math.PI * dof, 0.5 * D) /
            Math.pow(determinant, 0.5);
}
 

@Test
public void testMatrixDeterminant(){
    OpValidationSuite.ignoreFailing();  //Gradient check failing

    Nd4j.getRandom().setSeed(12345);
    INDArray in = Nd4j.rand(3,3);

    SameDiff sd = SameDiff.create();
    SDVariable var = sd.var("in", in);
    SDVariable md = sd.math().matrixDeterminant(var);

    double d = new LUDecomposition(CheckUtil.convertToApacheMatrix(in)).getDeterminant();


    INDArray outExp = Nd4j.scalar(d);

    String err = OpValidation.validate(new TestCase(sd)
            .expected(md.name(), outExp));
    assertNull(err);
}
 

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an symmetric n x n - matrix given as double[n][n]</li>
 * <li>b is an n - vector given as double[n],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 *
 * @param matrix The matrix A (left hand side of the linear equation).
 * @param vector The vector b (right hand of the linear equation).
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquationSymmetric(final double[][] matrix, final double[] vector) {
	if(isSolverUseApacheCommonsMath) {
		final DecompositionSolver solver = new LUDecomposition(new Array2DRowRealMatrix(matrix)).getSolver();
		return solver.solve(new Array2DRowRealMatrix(vector)).getColumn(0);
	}
	else {
		return org.jblas.Solve.solveSymmetric(new org.jblas.DoubleMatrix(matrix), new org.jblas.DoubleMatrix(vector)).data;
		/* To use the linear algebra package colt from cern.
		cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra();
		double[] x = linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray();

		return x;
		 */
	}
}
 

@Test
public void testInverseComparison() {

    List<Pair<INDArray, String>> list = NDArrayCreationUtil.getAllTestMatricesWithShape(10, 10, 12345, DataType.DOUBLE);

    for (Pair<INDArray, String> p : list) {
        INDArray orig = p.getFirst();
        orig.assign(Nd4j.rand(orig.shape()));
        INDArray inverse = InvertMatrix.invert(orig, false);
        RealMatrix rm = CheckUtil.convertToApacheMatrix(orig);
        RealMatrix rmInverse = new LUDecomposition(rm).getSolver().getInverse();

        INDArray expected = CheckUtil.convertFromApacheMatrix(rmInverse, orig.dataType());
        assertTrue(p.getSecond(), CheckUtil.checkEntries(expected, inverse, 1e-3, 1e-4));
    }
}
 

/**
 * Create a new East North Up system.
 * 
 * @param baseCoordinateLL the WGS84 coordinate to use a origin of the ENU. 
 */
public ENU( Coordinate baseCoordinateLL ) {
    checkZ(baseCoordinateLL);
    this._baseCoordinateLL = baseCoordinateLL;
    _lambda = toRadians(baseCoordinateLL.y);
    _phi = toRadians(baseCoordinateLL.x);
    _sinLambda = sin(_lambda);
    _cosLambda = cos(_lambda);
    _cosPhi = cos(_phi);
    _sinPhi = sin(_phi);
    _N = semiMajorAxis / sqrt(1 - eccentricityP2 * pow(_sinLambda, 2.0));

    double[][] rot = new double[][]{//
            {-_sinPhi, _cosPhi, 0}, //
            {-_cosPhi * _sinLambda, -_sinLambda * _sinPhi, _cosLambda}, //
            {_cosLambda * _cosPhi, _cosLambda * _sinPhi, _sinLambda},//
    };
    _rotationMatrix = MatrixUtils.createRealMatrix(rot);
    _inverseRotationMatrix = new LUDecomposition(_rotationMatrix).getSolver().getInverse();

    // the origin of the LTP expressed in ECEF-r coordinates
    double h = _baseCoordinateLL.z;
    _ecefROriginX = (h + _N) * _cosLambda * _cosPhi;
    _ecefROriginY = (h + _N) * _cosLambda * _sinPhi;
    _ecefROriginZ = (h + (1 - eccentricityP2) * _N) * _sinLambda;
}
 

@Test
public void testInverseComparison() {

    List<Pair<INDArray, String>> list = NDArrayCreationUtil.getAllTestMatricesWithShape(10, 10, 12345);

    for (Pair<INDArray, String> p : list) {
        INDArray orig = p.getFirst();
        orig.assign(Nd4j.rand(orig.shape()));
        INDArray inverse = InvertMatrix.invert(orig, false);
        RealMatrix rm = CheckUtil.convertToApacheMatrix(orig);
        RealMatrix rmInverse = new LUDecomposition(rm).getSolver().getInverse();

        INDArray expected = CheckUtil.convertFromApacheMatrix(rmInverse);
        assertTrue(p.getSecond(), CheckUtil.checkEntries(expected, inverse, 1e-3, 1e-4));
    }
}
 

/**
 * @see AbstractSolver#solve(AbstractMatrix, AbstractVector)
 */
@Override
public AbstractVector solve(AbstractMatrix m, AbstractVector b) {
  if (m instanceof ApacheMatrix && b instanceof ApacheVector) {
    DecompositionSolver solver = new LUDecomposition(((ApacheMatrix) m).getMatrix()).getSolver();
    RealVector bApache = ((ApacheVector) b).getVector();
    RealVector solution = solver.solve(bApache);
    return new ApacheVector(solution);
  }
  throw new UnsupportedOperationException();
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}
 

/**
 * Get the inverse of the covariance.
 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
 * @return inverse of the covariance
 */
protected RealMatrix getOmegaInverse() {
    if (OmegaInverse == null) {
        OmegaInverse = new LUDecomposition(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}
 

/**
 * Inverts a matrix
 * @param arr the array to invert
 * @param inPlace Whether to store the result in {@code arr}
 * @return the inverted matrix
 */
public static INDArray invert(INDArray arr, boolean inPlace) {
    if(arr.rank() == 2 && arr.length() == 1){
        //[1,1] edge case. Matrix inversion: [x] * [1/x] = [1]
        if(inPlace){
            return arr.rdivi(1.0);
        } else {
            return arr.rdiv(1.0);
        }
    }
    if (!arr.isSquare()) {
        throw new IllegalArgumentException("invalid array: must be square matrix");
    }

    //FIX ME: Please
   /* int[] IPIV = new int[arr.length() + 1];
    int LWORK = arr.length() * arr.length();
    INDArray WORK = Nd4j.create(new double[LWORK]);
    INDArray inverse = inPlace ? arr : arr.dup();
    Nd4j.getBlasWrapper().lapack().getrf(arr);
    Nd4j.getBlasWrapper().lapack().getri(arr.size(0),inverse,arr.size(0),IPIV,WORK,LWORK,0);*/

    RealMatrix rm = CheckUtil.convertToApacheMatrix(arr);
    RealMatrix rmInverse = new LUDecomposition(rm).getSolver().getInverse();


    INDArray inverse = CheckUtil.convertFromApacheMatrix(rmInverse, arr.dataType());
    if (inPlace)
        arr.assign(inverse);
    return inverse;

}
 

/**
 * Constructor
 * @param matrix    the matrix to decompose
 */
private XLUD(RealMatrix matrix) {
    this.lud = new LUDecomposition(matrix);
    this.l = LazyValue.of(() -> toDataFrame(lud.getL()));
    this.u = LazyValue.of(() -> toDataFrame(lud.getU()));
    this.p = LazyValue.of(() -> toDataFrame(lud.getP()));
}
 

/**
 * Solve a linear matrix equation, or system of linear scalar equations.
 *
 * @param a Coefficient matrix.
 * @param b Ordinate or “dependent variable” values.
 * @return Solution to the system a x = b. Returned shape is identical to b.
 */
public static Array solve(Array a, Array b) {
    Array r = Array.factory(DataType.DOUBLE, b.getShape());
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray_Double(a);
    RealMatrix coefficients = new Array2DRowRealMatrix(aa, false);
    DecompositionSolver solver = new LUDecomposition(coefficients).getSolver();
    double[] bb = (double[]) ArrayUtil.copyToNDJavaArray_Double(b);
    RealVector constants = new ArrayRealVector(bb, false);
    RealVector solution = solver.solve(constants);
    for (int i = 0; i < r.getSize(); i++) {
        r.setDouble(i, solution.getEntry(i));
    }

    return r;
}
 

/**
 * Get the inverse of the covariance.
 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
 * @return inverse of the covariance
 */
protected RealMatrix getOmegaInverse() {
    if (OmegaInverse == null) {
        OmegaInverse = new LUDecomposition(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}
 

/**
 * Get the inverse of the covariance.
 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
 * @return inverse of the covariance
 */
protected RealMatrix getOmegaInverse() {
    if (OmegaInverse == null) {
        OmegaInverse = new LUDecomposition(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}
 

/**
 * Creates an instance.
 * 
 * @param lu The result of the LU decomposition, not null. $\mathbf{L}$ cannot be singular.
 */
public LUDecompositionCommonsResult(LUDecomposition lu) {
  ArgChecker.notNull(lu, "LU decomposition");
  ArgChecker.notNull(lu.getL(), "Matrix is singular; could not perform LU decomposition");
  _determinant = lu.getDeterminant();
  _l = CommonsMathWrapper.unwrap(lu.getL());
  _p = CommonsMathWrapper.unwrap(lu.getP());
  _pivot = lu.getPivot();
  _solver = lu.getSolver();
  _u = CommonsMathWrapper.unwrap(lu.getU());
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}
 

/**
 * Get the inverse of the covariance.
 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
 * @return inverse of the covariance
 */
protected RealMatrix getOmegaInverse() {
    if (OmegaInverse == null) {
        OmegaInverse = new LUDecomposition(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}
 

/**
 * Computes model parameters and parameter variance using a QR decomposition of the X matrix
 * @param y     the response vector
 * @param x     the design matrix
 */
private RealMatrix computeBetaQR(RealVector y, RealMatrix x) {
    final int n = x.getRowDimension();
    final int p = x.getColumnDimension();
    final int offset = hasIntercept() ? 1 : 0;
    final QRDecomposition decomposition = new QRDecomposition(x, threshold);
    final RealVector betaVector = decomposition.getSolver().solve(y);
    final RealVector residuals = y.subtract(x.operate(betaVector));
    this.rss = residuals.dotProduct(residuals);
    this.errorVariance = rss / (n - p);
    this.stdError = Math.sqrt(errorVariance);
    this.residuals = createResidualsFrame(residuals);
    final RealMatrix rAug = decomposition.getR().getSubMatrix(0, p - 1, 0, p - 1);
    final RealMatrix rInv = new LUDecomposition(rAug).getSolver().getInverse();
    final RealMatrix covMatrix = rInv.multiply(rInv.transpose()).scalarMultiply(errorVariance);
    final RealMatrix result = new Array2DRowRealMatrix(p, 2);
    if (hasIntercept()) {
        result.setEntry(0, 0, betaVector.getEntry(0));      //Intercept coefficient
        result.setEntry(0, 1, covMatrix.getEntry(0, 0));    //Intercept variance
    }
    for (int i = 0; i < getRegressors().size(); i++) {
        final int index = i + offset;
        final double variance = covMatrix.getEntry(index, index);
        result.setEntry(index, 1, variance);
        result.setEntry(index, 0, betaVector.getEntry(index));
    }
    return result;
}
 

/**
 * Get the inverse of the covariance.
 * <p>The inverse of the covariance matrix is lazily evaluated and cached.</p>
 * @return inverse of the covariance
 */
protected RealMatrix getOmegaInverse() {
    if (OmegaInverse == null) {
        OmegaInverse = new LUDecomposition(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}
 

@Test(dataProvider = "styles")
public void testInverse(DataFrameAlgebra.Lib lib, boolean parallel) {
    DataFrameAlgebra.LIBRARY.set(lib);
    Array.of(20, 77, 95, 135, 233, 245).forEach(count -> {
        final DataFrame<Integer,Integer> frame = random(count, count, parallel, double.class);
        final DataFrame<Integer,Integer> inverse = frame.inverse();
        final RealMatrix matrix = new LUDecomposition(toMatrix(frame)).getSolver().getInverse();
        assertEquals(inverse, matrix);
    });
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}
 

/**
 * Calculates beta by GLS.
 * <pre>
 *  b=(X' Omega^-1 X)^-1X'Omega^-1 y
 * </pre>
 * @return beta
 */
@Override
protected RealVector calculateBeta() {
    RealMatrix OI = getOmegaInverse();
    RealMatrix XT = getX().transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(getX());
    RealMatrix inverse = new LUDecomposition(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(getY());
}
 

@Test
public void testMatrixDeterminant3(){
    OpValidationSuite.ignoreFailing();  //Gradient checks failing
    Nd4j.getRandom().setSeed(12345);
    INDArray in = Nd4j.rand(3,3);
    //System.out.println(in.shapeInfoToString());   //Rank: 2,Offset: 0 Order: c Shape: [3,3],  stride: [3,1]
    //System.out.println(Arrays.toString(in.data().asFloat())); //[0.27620894, 0.21801452, 0.062078513, 7.348895E-4, 0.24149609, 0.4948205, 0.93483436, 0.52035654, 0.30292067]

    SameDiff sd = SameDiff.create();
    SDVariable var = sd.var("in", in);
    SDVariable md = sd.math().matrixDeterminant(var);

    double d = new LUDecomposition(CheckUtil.convertToApacheMatrix(in)).getDeterminant();

    //https://en.wikipedia.org/wiki/Determinant
    double[][] a = in.toDoubleMatrix();
    double d2 = a[0][0] * a[1][1] * a[2][2]
            + a[0][1] * a[1][2] * a[2][0]
            + a[0][2] * a[1][0] * a[2][1]
            - a[0][2] * a[1][1] * a[2][0]
            - a[0][1] * a[1][0] * a[2][2]
            - a[0][0] * a[1][2] * a[2][1];
    assertEquals(d, d2, 1e-6);          //Manual calc and Apache commons both match:    0.03589524995561552

    INDArray outExp = Nd4j.scalar(d);

    String err = OpValidation.validate(new TestCase(sd)
            .expected(md.name(), outExp));
    assertNull(err);
}