org.apache.commons.math.linear.LUDecompositionImpl Java Examples

/**
  * http://www.slideshare.net/fullscreen/srowen/matrix-factorization/16 
  * @param recentitemInteractions
  * @param productFeaturesInverse
  * @param idMap
  * @return
  */
 private double[][] computeUserFoldInMatrix(double[][] itemFactors) 
 {
 	try
 	{
 		RealMatrix Y = new Array2DRowRealMatrix(itemFactors);
 		RealMatrix YTY = Y.transpose().multiply(Y);
 		RealMatrix YTYInverse = new LUDecompositionImpl(YTY).getSolver().getInverse();

 		return Y.multiply(YTYInverse).getData();
 	}
 	catch (InvalidMatrixException e)
 	{
 		logger.warn("Failed to create inverse of products feature matrix",e);
 		return null;
 	}
}
 

/**
 * Get the covariance matrix of the optimized parameters.
 *
 * @return the covariance matrix.
 * @throws org.apache.commons.math.linear.SingularMatrixException
 * if the covariance matrix cannot be computed (singular problem).
 * @throws org.apache.commons.math.exception.MathUserException if the
 * jacobian function throws one.
 */
public double[][] getCovariances() {
    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < rows; ++k) {
                sum += weightedResidualJacobian[k][i] * weightedResidualJacobian[k][j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    // compute the covariances matrix
    RealMatrix inverse =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
    return inverse.getData();
}
 

/**
 * Get the covariance matrix of unbound estimated parameters.
 * @param problem estimation problem
 * @return covariance matrix
 * @exception EstimationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances(EstimationProblem problem)
  throws EstimationException {

    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    final int n = problem.getMeasurements().length;
    final int m = problem.getUnboundParameters().length;
    final int max  = m * n;
    double[][] jTj = new double[m][m];
    for (int i = 0; i < m; ++i) {
        for (int j = i; j < m; ++j) {
            double sum = 0;
            for (int k = 0; k < max; k += m) {
                sum += jacobian[k + i] * jacobian[k + j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        RealMatrix inverse =
            new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
        return inverse.getData();
    } catch (InvalidMatrixException ime) {
        throw new EstimationException("unable to compute covariances: singular problem");
    }

}
 

/**
 * 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 = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}
 

/** test that U is upper triangular */
public void testUUpperTriangular() {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    RealMatrix u = new LUDecompositionImpl(matrix).getU();
    for (int i = 0; i < u.getRowDimension(); i++) {
        for (int j = 0; j < i; j++) {
            assertEquals(u.getEntry(i, j), 0, entryTolerance);
        }
    }
}
 

/**
 * 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 LUDecompositionImpl(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 LUDecompositionImpl(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 = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}
 

/** test matrices values */
public void testMatricesValues2() {
   LUDecomposition lu =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(luData));
    RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
            {    1.0,    0.0, 0.0 },
            {    0.0,    1.0, 0.0 },
            { 1.0 / 3.0, 0.0, 1.0 }
    });
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 6.0, 9.0,    8.0    },
            { 0.0, 5.0,    7.0    },
            { 0.0, 0.0, 1.0 / 3.0 }
    });
    RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
            { 0.0, 0.0, 1.0 },
            { 0.0, 1.0, 0.0 },
            { 1.0, 0.0, 0.0 }
    });
    int[] pivotRef = { 2, 1, 0 };

    // check values against known references
    RealMatrix l = lu.getL();
    assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
    RealMatrix u = lu.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
    RealMatrix p = lu.getP();
    assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
    int[] pivot = lu.getPivot();
    for (int i = 0; i < pivotRef.length; ++i) {
        assertEquals(pivotRef[i], pivot[i]);
    }

    // check the same cached instance is returned the second time
    assertTrue(l == lu.getL());
    assertTrue(u == lu.getU());
    assertTrue(p == lu.getP());
    
}
 

/**
 * 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 = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}
 

/**
 * 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 LUDecompositionImpl(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 = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}
 

/** test solve */
public void testSolve() {
    DecompositionSolver solver =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            { 1, 0 }, { 2, -5 }, { 3, 1 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { 19, -71 }, { -6, 22 }, { -2, 9 }
    });

    // using RealMatrix
    assertEquals(0, solver.solve(b).subtract(xRef).getNorm(), 1.0e-13);

    // using double[]
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     new ArrayRealVector(solver.solve(b.getColumn(i))).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using ArrayRealVector
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        assertEquals(0,
                     solver.solve(b.getColumnVector(i)).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

    // using RealVector with an alternate implementation
    for (int i = 0; i < b.getColumnDimension(); ++i) {
        ArrayRealVectorTest.RealVectorTestImpl v =
            new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(i));
        assertEquals(0,
                     solver.solve(v).subtract(xRef.getColumnVector(i)).getNorm(),
                     1.0e-13);
    }

}
 

/**
 * 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 = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}
 

/**
 * Get the covariance matrix of optimized parameters.
 * @return covariance matrix
 * @exception FunctionEvaluationException if the function jacobian cannot
 * be evaluated
 * @exception OptimizationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances()
    throws FunctionEvaluationException, OptimizationException {

    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < rows; ++k) {
                sum += jacobian[k][i] * jacobian[k][j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        RealMatrix inverse =
            new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
        return inverse.getData();
    } catch (InvalidMatrixException ime) {
        throw new OptimizationException("unable to compute covariances: singular problem");
    }

}
 

/**
 * 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 LUDecompositionImpl(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 LUDecompositionImpl(Omega).getSolver().getInverse();
    }
    return OmegaInverse;
}
 

/** test matrices values */
public void testMatricesValues1() {
   LUDecomposition lu =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData));
    RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0, 0.0, 0.0 },
            { 0.5, 1.0, 0.0 },
            { 0.5, 0.2, 1.0 }
    });
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 2.0,  5.0, 3.0 },
            { 0.0, -2.5, 6.5 },
            { 0.0,  0.0, 0.2 }
    });
    RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
            { 0.0, 1.0, 0.0 },
            { 0.0, 0.0, 1.0 },
            { 1.0, 0.0, 0.0 }
    });
    int[] pivotRef = { 1, 2, 0 };

    // check values against known references
    RealMatrix l = lu.getL();
    assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
    RealMatrix u = lu.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
    RealMatrix p = lu.getP();
    assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
    int[] pivot = lu.getPivot();
    for (int i = 0; i < pivotRef.length; ++i) {
        assertEquals(pivotRef[i], pivot[i]);
    }

    // check the same cached instance is returned the second time
    assertTrue(l == lu.getL());
    assertTrue(u == lu.getU());
    assertTrue(p == lu.getP());
    
}
 

/** test threshold impact */
public void testThreshold() {
    final RealMatrix matrix = MatrixUtils.createRealMatrix(new double[][] {
                                                   { 1.0, 2.0, 3.0},
                                                   { 2.0, 5.0, 3.0},
                                                   { 4.000001, 9.0, 9.0}
                                                 });
    assertFalse(new LUDecompositionImpl(matrix, 1.0e-5).getSolver().isNonSingular());
    assertTrue(new LUDecompositionImpl(matrix, 1.0e-10).getSolver().isNonSingular());
}
 

/** test that L is lower triangular with unit diagonal */
public void testLLowerTriangular() {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    RealMatrix l = new LUDecompositionImpl(matrix).getL();
    for (int i = 0; i < l.getRowDimension(); i++) {
        assertEquals(l.getEntry(i, i), 1, entryTolerance);
        for (int j = i + 1; j < l.getColumnDimension(); j++) {
            assertEquals(l.getEntry(i, j), 0, entryTolerance);
        }
    }
}
 

/** test that U is upper triangular */
public void testUUpperTriangular() {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    RealMatrix u = new LUDecompositionImpl(matrix).getU();
    for (int i = 0; i < u.getRowDimension(); i++) {
        for (int j = 0; j < i; j++) {
            assertEquals(u.getEntry(i, j), 0, entryTolerance);
        }
    }
}
 

/** test dimensions */
public void testDimensions() {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    LUDecomposition LU = new LUDecompositionImpl(matrix);
    assertEquals(testData.length, LU.getL().getRowDimension());
    assertEquals(testData.length, LU.getL().getColumnDimension());
    assertEquals(testData.length, LU.getU().getRowDimension());
    assertEquals(testData.length, LU.getU().getColumnDimension());
    assertEquals(testData.length, LU.getP().getRowDimension());
    assertEquals(testData.length, LU.getP().getColumnDimension());

}
 

/** test non-square matrix */
public void testNonSquare() {
    try {
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(new double[3][2]));
    } catch (InvalidMatrixException ime) {
        // expected behavior
    } catch (Exception e) {
        fail("wrong exception caught");
    }
}
 

/**
 * Get the covariance matrix of optimized parameters.
 * @return covariance matrix
 * @exception FunctionEvaluationException if the function jacobian cannot
 * be evaluated
 * @exception OptimizationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances()
    throws FunctionEvaluationException, OptimizationException {

    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    double[][] jTj = new double[cols][cols];
    for (int i = 0; i < cols; ++i) {
        for (int j = i; j < cols; ++j) {
            double sum = 0;
            for (int k = 0; k < rows; ++k) {
                sum += jacobian[k][i] * jacobian[k][j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        RealMatrix inverse =
            new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
        return inverse.getData();
    } catch (InvalidMatrixException ime) {
        throw new OptimizationException("unable to compute covariances: singular problem");
    }

}
 

/** test singular */
public void testSingular() {
    DecompositionSolver solver =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
    assertTrue(solver.isNonSingular());
    solver = new LUDecompositionImpl(MatrixUtils.createRealMatrix(singular)).getSolver();
    assertFalse(solver.isNonSingular());
    solver = new LUDecompositionImpl(MatrixUtils.createRealMatrix(bigSingular)).getSolver();
    assertFalse(solver.isNonSingular());
}
 

/** test matrices values */
public void testMatricesValues2() {
   LUDecomposition lu =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(luData));
    RealMatrix lRef = MatrixUtils.createRealMatrix(new double[][] {
            {    1.0,    0.0, 0.0 },
            {    0.0,    1.0, 0.0 },
            { 1.0 / 3.0, 0.0, 1.0 }
    });
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            { 6.0, 9.0,    8.0    },
            { 0.0, 5.0,    7.0    },
            { 0.0, 0.0, 1.0 / 3.0 }
    });
    RealMatrix pRef = MatrixUtils.createRealMatrix(new double[][] {
            { 0.0, 0.0, 1.0 },
            { 0.0, 1.0, 0.0 },
            { 1.0, 0.0, 0.0 }
    });
    int[] pivotRef = { 2, 1, 0 };

    // check values against known references
    RealMatrix l = lu.getL();
    assertEquals(0, l.subtract(lRef).getNorm(), 1.0e-13);
    RealMatrix u = lu.getU();
    assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-13);
    RealMatrix p = lu.getP();
    assertEquals(0, p.subtract(pRef).getNorm(), 1.0e-13);
    int[] pivot = lu.getPivot();
    for (int i = 0; i < pivotRef.length; ++i) {
        assertEquals(pivotRef[i], pivot[i]);
    }

    // check the same cached instance is returned the second time
    assertTrue(l == lu.getL());
    assertTrue(u == lu.getU());
    assertTrue(p == lu.getP());
    
}
 

/**
 * Get the covariance matrix of unbound estimated parameters.
 * @param problem estimation problem
 * @return covariance matrix
 * @exception EstimationException if the covariance matrix
 * cannot be computed (singular problem)
 */
public double[][] getCovariances(EstimationProblem problem)
  throws EstimationException {

    // set up the jacobian
    updateJacobian();

    // compute transpose(J).J, avoiding building big intermediate matrices
    final int n = problem.getMeasurements().length;
    final int m = problem.getUnboundParameters().length;
    final int max  = m * n;
    double[][] jTj = new double[m][m];
    for (int i = 0; i < m; ++i) {
        for (int j = i; j < m; ++j) {
            double sum = 0;
            for (int k = 0; k < max; k += m) {
                sum += jacobian[k + i] * jacobian[k + j];
            }
            jTj[i][j] = sum;
            jTj[j][i] = sum;
        }
    }

    try {
        // compute the covariances matrix
        RealMatrix inverse =
            new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
        return inverse.getData();
    } catch (InvalidMatrixException ime) {
        throw new EstimationException("unable to compute covariances: singular problem");
    }

}
 

/** test singular */
public void testSingular() {
    DecompositionSolver solver =
        new LUDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
    assertTrue(solver.isNonSingular());
    solver = new LUDecompositionImpl(MatrixUtils.createRealMatrix(singular)).getSolver();
    assertFalse(solver.isNonSingular());
    solver = new LUDecompositionImpl(MatrixUtils.createRealMatrix(bigSingular)).getSolver();
    assertFalse(solver.isNonSingular());
}
 

/**
 * 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 LUDecompositionImpl(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 = X.transpose();
    RealMatrix XTOIX = XT.multiply(OI).multiply(X);
    RealMatrix inverse = new LUDecompositionImpl(XTOIX).getSolver().getInverse();
    return inverse.multiply(XT).multiply(OI).operate(Y);
}