Java Code Examples for org.apache.commons.math.linear.RealMatrix#getColumnDimension()

/**
 * <p>Compute the "hat" matrix.
 * </p>
 * <p>The hat matrix is defined in terms of the design matrix X
 *  by X(X<sup>T</sup>X)<sup>-1</sup>X<sup>T</sup>
 * </p>
 * <p>The implementation here uses the QR decomposition to compute the
 * hat matrix as Q I<sub>p</sub>Q<sup>T</sup> where I<sub>p</sub> is the
 * p-dimensional identity matrix augmented by 0's.  This computational
 * formula is from "The Hat Matrix in Regression and ANOVA",
 * David C. Hoaglin and Roy E. Welsch,
 * <i>The American Statistician</i>, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
 *
 * @return the hat matrix
 */
public RealMatrix calculateHat() {
    // Create augmented identity matrix
    RealMatrix Q = qr.getQ();
    final int p = qr.getR().getColumnDimension();
    final int n = Q.getColumnDimension();
    Array2DRowRealMatrix augI = new Array2DRowRealMatrix(n, n);
    double[][] augIData = augI.getDataRef();
    for (int i = 0; i < n; i++) {
        for (int j =0; j < n; j++) {
            if (i == j && i < p) {
                augIData[i][j] = 1d;
            } else {
                augIData[i][j] = 0d;
            }
        }
    }

    // Compute and return Hat matrix
    return Q.multiply(augI).multiply(Q.transpose());
}
 

/**
 * <p>Compute the "hat" matrix.
 * </p>
 * <p>The hat matrix is defined in terms of the design matrix X
 *  by X(X<sup>T</sup>X)<sup>-1</sup>X<sup>T</sup>
 * </p>
 * <p>The implementation here uses the QR decomposition to compute the
 * hat matrix as Q I<sub>p</sub>Q<sup>T</sup> where I<sub>p</sub> is the
 * p-dimensional identity matrix augmented by 0's.  This computational
 * formula is from "The Hat Matrix in Regression and ANOVA",
 * David C. Hoaglin and Roy E. Welsch,
 * <i>The American Statistician</i>, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
 *
 * @return the hat matrix
 */
public RealMatrix calculateHat() {
    // Create augmented identity matrix
    RealMatrix Q = qr.getQ();
    final int p = qr.getR().getColumnDimension();
    final int n = Q.getColumnDimension();
    Array2DRowRealMatrix augI = new Array2DRowRealMatrix(n, n);
    double[][] augIData = augI.getDataRef();
    for (int i = 0; i < n; i++) {
        for (int j =0; j < n; j++) {
            if (i == j && i < p) {
                augIData[i][j] = 1d;
            } else {
                augIData[i][j] = 0d;
            }
        }
    }

    // Compute and return Hat matrix
    return Q.multiply(augI).multiply(Q.transpose());
}
 

/**
 * Derives a correlation matrix from a covariance matrix.
 *
 * <p>Uses the formula <br/>
 * <code>r(X,Y) = cov(X,Y)/s(X)s(Y)</code> where
 * <code>r(&middot,&middot;)</code> is the correlation coefficient and
 * <code>s(&middot;)</code> means standard deviation.</p>
 *
 * @param covarianceMatrix the covariance matrix
 * @return correlation matrix
 */
public RealMatrix covarianceToCorrelation(RealMatrix covarianceMatrix) {
    int nVars = covarianceMatrix.getColumnDimension();
    RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
    for (int i = 0; i < nVars; i++) {
        double sigma = Math.sqrt(covarianceMatrix.getEntry(i, i));
        outMatrix.setEntry(i, i, 1d);
        for (int j = 0; j < i; j++) {
            double entry = covarianceMatrix.getEntry(i, j) /
                   (sigma * Math.sqrt(covarianceMatrix.getEntry(j, j)));
            outMatrix.setEntry(i, j, entry);
            outMatrix.setEntry(j, i, entry);
        }
    }
    return outMatrix;
}
 

/**
 * Returns true iff there is a column that is a scalar multiple of column
 * in searchMatrix (modulo tolerance)
 */
private boolean isIncludedColumn(double[] column, RealMatrix searchMatrix,
        double tolerance) {
    boolean found = false;
    int i = 0;
    while (!found && i < searchMatrix.getColumnDimension()) {
        double multiplier = 1.0;
        boolean matching = true;
        int j = 0;
        while (matching && j < searchMatrix.getRowDimension()) {
            double colEntry = searchMatrix.getEntry(j, i);
            // Use the first entry where both are non-zero as scalar
            if (Math.abs(multiplier - 1.0) <= Math.ulp(1.0) && Math.abs(colEntry) > 1E-14
                    && Math.abs(column[j]) > 1e-14) {
                multiplier = colEntry / column[j];
            }
            if (Math.abs(column[j] * multiplier - colEntry) > tolerance) {
                matching = false;
            }
            j++;
        }
        found = matching;
        i++;
    }
    return found;
}
 

/**
 * Derives a correlation matrix from a covariance matrix.
 *
 * <p>Uses the formula <br/>
 * <code>r(X,Y) = cov(X,Y)/s(X)s(Y)</code> where
 * <code>r(&middot,&middot;)</code> is the correlation coefficient and
 * <code>s(&middot;)</code> means standard deviation.</p>
 *
 * @param covarianceMatrix the covariance matrix
 * @return correlation matrix
 */
public RealMatrix covarianceToCorrelation(RealMatrix covarianceMatrix) {
    int nVars = covarianceMatrix.getColumnDimension();
    RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
    for (int i = 0; i < nVars; i++) {
        double sigma = Math.sqrt(covarianceMatrix.getEntry(i, i));
        outMatrix.setEntry(i, i, 1d);
        for (int j = 0; j < i; j++) {
            double entry = covarianceMatrix.getEntry(i, j) /
                   (sigma * Math.sqrt(covarianceMatrix.getEntry(j, j)));
            outMatrix.setEntry(i, j, entry);
            outMatrix.setEntry(j, i, entry);
        }
    }
    return outMatrix;
}
 

private void checkTriDiagonal(RealMatrix m) {
    final int rows = m.getRowDimension();
    final int cols = m.getColumnDimension();
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            if ((i < j - 1) || (i > j + 1)) {
                assertEquals(0, m.getEntry(i, j), 1.0e-16);
            }
        }
    }
}
 

/**
 * Compute a covariance matrix from a matrix whose columns represent
 * covariates.
 * @param matrix input matrix (must have at least two columns and two rows)
 * @param biasCorrected determines whether or not covariance estimates are bias-corrected
 * @return covariance matrix
 */
protected RealMatrix computeCovarianceMatrix(RealMatrix matrix, boolean biasCorrected) {
    int dimension = matrix.getColumnDimension();
    Variance variance = new Variance(biasCorrected);
    RealMatrix outMatrix = new BlockRealMatrix(dimension, dimension);
    for (int i = 0; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
          double cov = covariance(matrix.getColumn(i), matrix.getColumn(j), biasCorrected);
          outMatrix.setEntry(i, j, cov);
          outMatrix.setEntry(j, i, cov);
        }
        outMatrix.setEntry(i, i, variance.evaluate(matrix.getColumn(i)));
    }
    return outMatrix;
}
 

private void checkDimension(RealMatrix m) {
    int rows = m.getRowDimension();
    int columns = m.getColumnDimension();
    QRDecomposition qr = new QRDecompositionImpl(m);
    assertEquals(rows,    qr.getQ().getRowDimension());
    assertEquals(rows,    qr.getQ().getColumnDimension());
    assertEquals(rows,    qr.getR().getRowDimension());
    assertEquals(columns, qr.getR().getColumnDimension());        
}
 

private void checkDimension(RealMatrix m) {
    int rows = m.getRowDimension();
    int columns = m.getColumnDimension();
    QRDecomposition qr = new QRDecompositionImpl(m);
    assertEquals(rows,    qr.getQ().getRowDimension());
    assertEquals(rows,    qr.getQ().getColumnDimension());
    assertEquals(rows,    qr.getR().getRowDimension());
    assertEquals(columns, qr.getR().getColumnDimension());
}
 

/**
 * Computes the correlation matrix for the columns of the
 * input matrix.
 *
 * @param matrix matrix with columns representing variables to correlate
 * @return correlation matrix
 */
public RealMatrix computeCorrelationMatrix(RealMatrix matrix) {
    int nVars = matrix.getColumnDimension();
    RealMatrix outMatrix = new BlockRealMatrix(nVars, nVars);
    for (int i = 0; i < nVars; i++) {
        for (int j = 0; j < i; j++) {
          double corr = correlation(matrix.getColumn(i), matrix.getColumn(j));
          outMatrix.setEntry(i, j, corr);
          outMatrix.setEntry(j, i, corr);
        }
        outMatrix.setEntry(i, i, 1d);
    }
    return outMatrix;
}
 

/**
 * @param m Input matrix.
 * @param k diagonal position.
 * @return Upper triangular part of matrix.
 */
private static RealMatrix triu(final RealMatrix m, int k) {
    double[][] d = new double[m.getRowDimension()][m.getColumnDimension()];
    for (int r = 0; r < m.getRowDimension(); r++) {
        for (int c = 0; c < m.getColumnDimension(); c++) {
            d[r][c] = r <= c - k ? m.getEntry(r, c) : 0;
        }
    }
    return new Array2DRowRealMatrix(d, false);
}
 

public CorrelatedRandomVectorGeneratorTest() {
    mean = new double[] { 0.0, 1.0, -3.0, 2.3 };

    RealMatrix b = MatrixUtils.createRealMatrix(4, 3);
    int counter = 0;
    for (int i = 0; i < b.getRowDimension(); ++i) {
        for (int j = 0; j < b.getColumnDimension(); ++j) {
            b.setEntry(i, j, 1.0 + 0.1 * ++counter);
        }
    }
    RealMatrix bbt = b.multiply(b.transpose());
    covariance = MatrixUtils.createRealMatrix(mean.length, mean.length);
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        covariance.setEntry(i, i, bbt.getEntry(i, i));
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            double s = bbt.getEntry(i, j);
            covariance.setEntry(i, j, s);
            covariance.setEntry(j, i, s);
        }
    }

    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(17399225432l);
    GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
    generator = new CorrelatedRandomVectorGenerator(mean,
                                                    covariance,
                                                    1.0e-12 * covariance.getNorm(),
                                                    rawGenerator);
}
 

/** test that L is lower triangular */
@Test
public void testLLowerTriangular() throws MathException {
    RealMatrix matrix = MatrixUtils.createRealMatrix(testData);
    RealMatrix l = new CholeskyDecompositionImpl(matrix).getL();
    for (int i = 0; i < l.getRowDimension(); i++) {
        for (int j = i + 1; j < l.getColumnDimension(); j++) {
            assertEquals(0.0, l.getEntry(i, j), 0.0);
        }
    }
}
 

/**
 * Compute a covariance matrix from a matrix whose columns represent
 * covariates.
 * @param matrix input matrix (must have at least two columns and two rows)
 * @param biasCorrected determines whether or not covariance estimates are bias-corrected
 * @return covariance matrix
 */
protected RealMatrix computeCovarianceMatrix(RealMatrix matrix, boolean biasCorrected) {
    int dimension = matrix.getColumnDimension();
    Variance variance = new Variance(biasCorrected);
    RealMatrix outMatrix = new BlockRealMatrix(dimension, dimension);
    for (int i = 0; i < dimension; i++) {
        for (int j = 0; j < i; j++) {
          double cov = covariance(matrix.getColumn(i), matrix.getColumn(j), biasCorrected);
          outMatrix.setEntry(i, j, cov);
          outMatrix.setEntry(j, i, cov);
        }
        outMatrix.setEntry(i, i, variance.evaluate(matrix.getColumn(i)));
    }
    return outMatrix;
}
 

/**
 * Throws IllegalArgumentException of the matrix does not have at least
 * two columns and two rows
 * @param matrix matrix to check
 */
private void checkSufficientData(final RealMatrix matrix) {
    int nRows = matrix.getRowDimension();
    int nCols = matrix.getColumnDimension();
    if (nRows < 2 || nCols < 2) {
        throw MathRuntimeException.createIllegalArgumentException(
                LocalizedFormats.INSUFFICIENT_ROWS_AND_COLUMNS,
                nRows, nCols);
    }
}
 

protected void fillUpper(RealMatrix matrix, double diagonalValue) {
    int dimension = matrix.getColumnDimension();
    for (int i = 0; i < dimension; i++) {
        matrix.setEntry(i, i, diagonalValue);
        for (int j = i+1; j < dimension; j++) {
            matrix.setEntry(i, j, matrix.getEntry(j, i));
        }
    }
}
 

public static RealMatrix expandWithZeros (RealMatrix original, int targetNumRows, int targetNumCols){
    int r, c;
    RealMatrix newMatrix = MatrixTools.cleanArray2DMatrix(targetNumRows, targetNumCols);
    for (r = 0; r < original.getRowDimension(); r++){
        for (c = 0; c < original.getColumnDimension(); c++){
            newMatrix.addToEntry(r, c, original.getEntry(r, c));
        }
    }
    return newMatrix;
}
 

/**
 * @param mat
 *            Input matrix.
 * @param n
 *            Number of row replicates.
 * @param m
 *            Number of column replicates.
 * @return Matrix which replicates the input matrix in both directions.
 */
private static RealMatrix repmat(final RealMatrix mat, int n, int m) {
    int rd = mat.getRowDimension();
    int cd = mat.getColumnDimension();
    double[][] d = new double[n * rd][m * cd];
    for (int r = 0; r < n * rd; r++) {
        for (int c = 0; c < m * cd; c++) {
            d[r][c] = mat.getEntry(r % rd, c % cd);
        }
    }
    return new Array2DRowRealMatrix(d, false);
}
 

/** test solve */
public void testSolve() throws MathException {
    DecompositionSolver solver =
        new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(testData)).getSolver();
    RealMatrix b = MatrixUtils.createRealMatrix(new double[][] {
            {   78,  -13,    1 },
            {  414,  -62,   -1 },
            { 1312, -202,  -37 },
            { 2989, -542,  145 },
            { 5510, -1465, 201 }
    });
    RealMatrix xRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1,  0,  1 },
            { 0,  1,  1 },
            { 2,  1, -4 },
            { 2,  2,  2 },
            { 5, -3,  0 }
    });

    // 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);
    }

}
 

/** Build a simple converter for correlated residuals with the specific weights.
 * <p>
 * The scalar objective function value is computed as:
 * <pre>
 * objective = y<sup>T</sup>y with y = scale&times;(observation-objective)
 * </pre>
 * </p>
 * <p>
 * The array computed by the objective function, the observations array and the
 * the scaling matrix must have consistent sizes or a {@link FunctionEvaluationException}
 * will be triggered while computing the scalar objective.
 * </p>
 * @param function vectorial residuals function to wrap
 * @param observations observations to be compared to objective function to compute residuals
 * @param scale scaling matrix
 * @exception IllegalArgumentException if the observations vector and the scale
 * matrix dimensions don't match (objective function dimension is checked only when
 * the {@link #value(double[])} method is called)
 */
public LeastSquaresConverter(final MultivariateVectorialFunction function,
                             final double[] observations, final RealMatrix scale)
    throws IllegalArgumentException {
    if (observations.length != scale.getColumnDimension()) {
        throw MathRuntimeException.createIllegalArgumentException(
                "dimension mismatch {0} != {1}",
                observations.length, scale.getColumnDimension());
    }
    this.function     = function;
    this.observations = observations.clone();
    this.weights      = null;
    this.scale        = scale.copy();
}