/** * Solve A*X = b * * @param b right hand side * @return solution if A is square, least squares solution otherwise */ public static double[] solve(double[][] A, double[] b) { final int rows = A.length, cols = A[0].length; return rows == cols // ? (new LUDecomposition(A, rows, cols)).solve(b) // : (new QRDecomposition(A, rows, cols)).solve(b); }
/** * Matrix determinant * * @return determinant */ public final double det() { return new LUDecomposition(this).det(); }
/** * Matrix inverse or pseudoinverse * * @param A matrix to invert * @return inverse(A) if A is square, pseudoinverse otherwise. */ public static double[][] inverse(double[][] A) { final int rows = A.length, cols = A[0].length; return rows == cols // ? (new LUDecomposition(A, rows, cols)).inverse() // : (new QRDecomposition(A, rows, cols)).inverse(); }
/** * Matrix inverse for square matrixes. * * @return inverse(A), or inverse(A + epsilon E) if singular. */ public final Matrix robustInverse() { LUDecomposition d = new LUDecomposition(this); if(!d.isNonsingular()) { d = new LUDecomposition(plusDiagonal(SINGULARITY_CHEAT).getArrayRef(), elements.length, columndimension); } return d.solve(identity(elements.length, elements.length)); }
/** * Computes the loglikelihood of all normal objects. Gaussian model * * @param objids Object IDs for 'normal' objects. * @param builder Covariance matrix builder * @param relation Database * @return loglikelihood for normal objects */ private double loglikelihoodNormal(DBIDs objids, SetDBIDs anomalous, CovarianceMatrix builder, Relation<V> relation) { double[] mean = builder.getMeanVector(); final LUDecomposition lu = new LUDecomposition(builder.makeSampleMatrix()); double[][] covInv = lu.inverse(); // for each object compute probability and sum double prob = (objids.size() - anomalous.size()) * -FastMath.log(FastMath.sqrt(MathUtil.powi(MathUtil.TWOPI, RelationUtil.dimensionality(relation)) * lu.det())); for(DBIDIter iter = objids.iter(); iter.valid(); iter.advance()) { if(!anomalous.contains(iter)) { double[] xcent = minusEquals(relation.get(iter).toArray(), mean); prob -= .5 * transposeTimesTimes(xcent, covInv, xcent); } } return prob; }
@Override public void finalizeEStep() { final int dim = mean.getDimensionality(); // TODO: improve handling of degenerated cases? if(wsum > Double.MIN_NORMAL) { covariance.timesEquals(1. / wsum); } LUDecomposition lu = new LUDecomposition(covariance); double det = lu.det(); if(!(det > 0.)) { // Add a small value to the diagonal covariance.plusDiagonalEquals(Matrix.SINGULARITY_CHEAT); lu = new LUDecomposition(covariance); // Should no longer be zero now. det = lu.det(); if(!(det > 0.)) { LOG.warning("Singularity cheat did not resolve zero determinant."); // assert (det > 0) : "Singularity cheat did not resolve zero determinant."; det = 1.; } } normDistrFactor = 1. / Math.sqrt(norm * det); invCovMatr = lu.solve(Matrix.identity(dim, dim)); }
@Test public void testWikipediaQR() { double[][] M = { // { 12, -51, 4 }, // { 6, 167, -68 }, // { -4, 24, -41 } }; LUDecomposition lu = new LUDecomposition(M); final double[][] l = lu.getL(), u = lu.getU(); checkTriangular(l, u); final double[][] B = times(l, u); // Rearrange final double[][] M2 = getMatrix(M, lu.getPivot(), 0, M.length); assertTrue("Not a proper decomposition.", almostEquals(M2, B, 1e-15)); }
/** * Testing the isNonSigular methods of the LUDecomposition class. */ @Test public void testIsNonSigular() { assertTrue(new LUDecomposition(M).isNonsingular()); assertFalse(new LUDecomposition(new double[][] { { 1, 1 }, { 1, 1 }, { 2, 2 } }).isNonsingular()); }
/** * Solve A*X = b * * @param b A column vector with as many rows as A * @return X so that L*U*X = b(piv) * @throws IllegalArgumentException Matrix row dimensions must agree. * @throws ArithmeticException Matrix is singular. */ public double[] solve(double[] b) { if(b.length != m) { throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS); } if(!this.isNonsingular()) { throw new ArithmeticException(ERR_SINGULAR); } double[] bc = new double[piv.length]; for(int i = 0; i < piv.length; i++) { bc[i] = b[piv[i]]; } solveInplace(bc); return n < bc.length ? Arrays.copyOf(bc, n) : bc; }
/** * Testing the getPivot and getDoublePivot methods of the LUDecomposition * class with {@link M} as data for testing. */ @Test public void testPivot() { // assert that octave and scipy function have almost same results assertTrue(almostEquals(U_OCTAVE, U_SCIPY, 1e-8)); // transformation of P_scipy into the format of the piv vector final int[] piv = { 7, 6, 9, 1, 10, 0, 4, 5, 8, 2, 3 }; assertArrayEquals(piv, new LUDecomposition(M).getPivot()); }
/** * Testing getU method of the LUDecomposition class with {@link M} as data for * testing. */ @Test public void testGetU() { // assert that octave and scipy function have almost same results assertTrue(almostEquals(U_OCTAVE, U_SCIPY, 1e-8)); final double[][] u = new LUDecomposition(M).getU(); assertTrue(almostEquals(u, U_OCTAVE, 1e-13)); assertTrue(almostEquals(u, U_SCIPY, 1e-8)); }
/** * Testing getL method of the LUDecomposition class with {@link M} as data for * testing. */ @Test public void testGetL() { // octave returns L in different format so we need to convert with times and // P_SCIPY assertTrue(almostEquals(L_OCTAVE, times(P_SCIPY, L_SCIPY), 1e-8)); double[][] m = new LUDecomposition(M).getL(); // assert that octave and scipy function have almost same results // octave returns L in different format so we need to convert with times and // P_SCIPY assertTrue(almostEquals(times(P_SCIPY, m), L_OCTAVE, 1e-9)); assertTrue(almostEquals(m, L_SCIPY, 1e-8)); }
/** * Solve A*X = b * * @param b A vector * @return b * @throws IllegalArgumentException Matrix row dimensions must agree. * @throws ArithmeticException Matrix is singular. */ public double[] solveInplace(double[] b) { if(b.length != m) { throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS); } if(!this.isNonsingular()) { throw new ArithmeticException(ERR_SINGULAR); } // Solve L*Y = B(piv,:) for(int k = 0; k < n; k++) { for(int i = k + 1; i < n; i++) { b[i] -= b[k] * LU[i][k]; } } // Solve U*X = Y; for(int k = n - 1; k >= 0; k--) { b[k] /= LU[k][k]; for(int i = 0; i < k; i++) { b[i] -= b[k] * LU[i][k]; } } return b; }
/** * Find the inverse matrix. * * @return Inverse matrix * @throws ArithmeticException Matrix is rank deficient. */ public double[][] inverse() { // Build permuted identity matrix efficiently: double[][] b = new double[piv.length][m]; for(int i = 0; i < piv.length; i++) { b[piv[i]][i] = 1.; } return solveInplace(b); } }
/** * Computes the loglikelihood of all normal objects. Gaussian model * * @param objids Object IDs for 'normal' objects. * @param builder Covariance matrix builder * @param relation Database * @return loglikelihood for normal objects */ private double loglikelihoodNormal(DBIDs objids, SetDBIDs anomalous, CovarianceMatrix builder, Relation<V> relation) { double[] mean = builder.getMeanVector(); final LUDecomposition lu = new LUDecomposition(builder.makeSampleMatrix()); double[][] covInv = lu.inverse(); // for each object compute probability and sum double prob = (objids.size() - anomalous.size()) * -FastMath.log(FastMath.sqrt(MathUtil.powi(MathUtil.TWOPI, RelationUtil.dimensionality(relation)) * lu.det())); for(DBIDIter iter = objids.iter(); iter.valid(); iter.advance()) { if(!anomalous.contains(iter)) { double[] xcent = minusEquals(relation.get(iter).toArray(), mean); prob -= .5 * transposeTimesTimes(xcent, covInv, xcent); } } return prob; }
@Test public void testJamaExample() { double[][] M = transpose(new double[][] { { 0., 2., 3. }, { 5., 6., 7. }, { 9., 10., 11. } }); LUDecomposition lu = new LUDecomposition(M); final double[][] l = lu.getL(), u = lu.getU(); checkTriangular(l, u); final double[][] B = times(l, u); // Rearrange final double[][] M2 = getMatrix(M, lu.getPivot(), 0, M.length); assertTrue("Not a proper decomposition.", almostEquals(M2, B, 1e-15)); }
/** * Solve A*X = b * * @param b A column vector with as many rows as A * @return X so that L*U*X = b(piv) * @throws IllegalArgumentException Matrix row dimensions must agree. * @throws ArithmeticException Matrix is singular. */ public double[] solve(double[] b) { if(b.length != m) { throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS); } if(!this.isNonsingular()) { throw new ArithmeticException(ERR_SINGULAR); } double[] bc = new double[piv.length]; for(int i = 0; i < piv.length; i++) { bc[i] = b[piv[i]]; } solveInplace(bc); return n < bc.length ? Arrays.copyOf(bc, n) : bc; }
/** * Solve A*X = b * * @param b A vector * @return b * @throws IllegalArgumentException Matrix row dimensions must agree. * @throws ArithmeticException Matrix is singular. */ public double[] solveInplace(double[] b) { if(b.length != m) { throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS); } if(!this.isNonsingular()) { throw new ArithmeticException(ERR_SINGULAR); } // Solve L*Y = B(piv,:) for(int k = 0; k < n; k++) { for(int i = k + 1; i < n; i++) { b[i] -= b[k] * LU[i][k]; } } // Solve U*X = Y; for(int k = n - 1; k >= 0; k--) { b[k] /= LU[k][k]; for(int i = 0; i < k; i++) { b[i] -= b[k] * LU[i][k]; } } return b; }
/** * Find the inverse matrix. * * @return Inverse matrix * @throws ArithmeticException Matrix is rank deficient. */ public double[][] inverse() { // Build permuted identity matrix efficiently: double[][] b = new double[piv.length][m]; for(int i = 0; i < piv.length; i++) { b[piv[i]][i] = 1.; } return solveInplace(b); } }
/** * Solve A*X = B * * @param B right hand side * @return solution if A is square, least squares solution otherwise */ public final Matrix solve(final Matrix B) { return (elements.length == columndimension ? (new LUDecomposition(this)).solve(B) : (new QRDecomposition(this)).solve(B)); }