Java Code Examples for org.apache.commons.math3.linear.RealMatrix#getSubMatrix()

/**
 * QR decomposition for a tridiagonal matrix T.
 *
 * @see https://gist.github.com/lightcatcher/8118181
 * @see http://www.ericmart.in/blog/optimizing_julia_tridiag_qr
 * @param T target tridiagonal matrix
 * @param R output matrix for R which is the same shape as T
 * @param Qt output matrix for Q.T which is the same shape an T
 */
public static void tridiagonalQR(@Nonnull final RealMatrix T, @Nonnull final RealMatrix R,
        @Nonnull final RealMatrix Qt) {
    int n = T.getRowDimension();
    Preconditions.checkArgument(n == R.getRowDimension() && n == R.getColumnDimension(),
        "T and R must be the same shape");
    Preconditions.checkArgument(n == Qt.getRowDimension() && n == Qt.getColumnDimension(),
        "T and Qt must be the same shape");

    // initial R = T
    R.setSubMatrix(T.getData(), 0, 0);

    // initial Qt = identity
    Qt.setSubMatrix(eye(n), 0, 0);

    for (int i = 0; i < n - 1; i++) {
        // Householder projection for a vector x
        // https://en.wikipedia.org/wiki/Householder_transformation
        RealVector x = T.getSubMatrix(i, i + 1, i, i).getColumnVector(0);
        x = unitL2norm(x);

        RealMatrix subR = R.getSubMatrix(i, i + 1, 0, n - 1);
        R.setSubMatrix(
            subR.subtract(x.outerProduct(subR.preMultiply(x)).scalarMultiply(2)).getData(), i,
            0);

        RealMatrix subQt = Qt.getSubMatrix(i, i + 1, 0, n - 1);
        Qt.setSubMatrix(
            subQt.subtract(x.outerProduct(subQt.preMultiply(x)).scalarMultiply(2)).getData(), i,
            0);
    }
}
 

/**
 * Asserts that a collection of beta-hats corresponds to the expected value given
 * the input pre-tangent normalization matrix and the PoN file.
 */
private void assertBetaHats(final ReadCountCollection preTangentNormalized,
                            final RealMatrix actual, final File ponFile) {
    Assert.assertEquals(actual.getColumnDimension(), preTangentNormalized.columnNames().size());
    final double epsilon = PCATangentNormalizationUtils.EPSILON;

    try (final HDF5File ponReader = new HDF5File(ponFile)) {
        final PCACoveragePoN pon = new HDF5PCACoveragePoN(ponReader);
        final List<String> ponTargets = pon.getPanelTargetNames();
        final RealMatrix inCounts = reorderTargetsToPoNOrder(preTangentNormalized, ponTargets);

        // obtain subset of relevant targets to calculate the beta-hats;
        final int[][] betaHatTargets = new int[inCounts.getColumnDimension()][];
        for (int i = 0; i < inCounts.getColumnDimension(); i++) {
            final List<Integer> relevantTargets = new ArrayList<>();
            for (int j = 0; j < inCounts.getRowDimension(); j++) {
                if (inCounts.getEntry(j, i) > 1  +  (Math.log(epsilon) / Math.log(2))) {
                    relevantTargets.add(j);
                }
            }
            betaHatTargets[i] = relevantTargets.stream().mapToInt(Integer::intValue).toArray();
        }
        // calculate beta-hats per column and check with actual values.
        final RealMatrix normalsInv = pon.getReducedPanelPInverseCounts();
        Assert.assertEquals(actual.getRowDimension(), normalsInv.getRowDimension());
        final RealMatrix normalsInvT = normalsInv.transpose();
        for (int i = 0; i < inCounts.getColumnDimension(); i++) {
            final RealMatrix inValues = inCounts.getColumnMatrix(i).transpose().getSubMatrix(new int[] { 0 }, betaHatTargets[i]);
            final RealMatrix normalValues = normalsInvT.getSubMatrix(betaHatTargets[i], IntStream.range(0, normalsInvT.getColumnDimension()).toArray());
            final RealMatrix betaHats = inValues.multiply(normalValues);
            for (int j = 0; j < actual.getRowDimension(); j++) {
                Assert.assertEquals(actual.getEntry(j, i), betaHats.getEntry(0, j),0.000001,"Col " + i + " row " + j);
            }
        }
    }
}
 

/**
 * @param x series of input in LIFO order
 * @param k AR window size
 * @return x_hat predicted x
 * @link https://en.wikipedia.org/wiki/Matrix_multiplication#Outer_product
 */
@Nonnull
public RealVector update(@Nonnull final ArrayRealVector[] x, final int k) {
    Preconditions.checkArgument(x.length >= 1, "x.length MUST be greater than 1: " + x.length);
    Preconditions.checkArgument(k >= 0, "k MUST be greater than or equals to 0: ", k);
    Preconditions.checkArgument(k < _C.length,
        "k MUST be less than |C| but " + "k=" + k + ", |C|=" + _C.length);

    final ArrayRealVector x_t = x[0];
    final int dims = x_t.getDimension();

    if (_initialized == false) {
        this._mu = x_t.copy();
        this._sigma = new BlockRealMatrix(dims, dims);
        assert (_sigma.isSquare());
        this._initialized = true;
        return new ArrayRealVector(dims);
    }
    Preconditions.checkArgument(k >= 1, "k MUST be greater than 0: ", k);

    // old parameters are accessible to compute the Hellinger distance
    this._muOld = _mu.copy();
    this._sigmaOld = _sigma.copy();

    // update mean vector
    // \hat{mu} := (1-r) \hat{µ} + r x_t
    this._mu = _mu.mapMultiply(1.d - _r).add(x_t.mapMultiply(_r));

    // compute residuals (x - \hat{µ})
    final RealVector[] xResidual = new RealVector[k + 1];
    for (int j = 0; j <= k; j++) {
        xResidual[j] = x[j].subtract(_mu);
    }

    // update covariance matrices
    // C_j := (1-r) C_j + r (x_t - \hat{µ}) (x_{t-j} - \hat{µ})'
    final RealMatrix[] C = this._C;
    final RealVector rxResidual0 = xResidual[0].mapMultiply(_r); // r (x_t - \hat{µ})
    for (int j = 0; j <= k; j++) {
        RealMatrix Cj = C[j];
        if (Cj == null) {
            C[j] = rxResidual0.outerProduct(x[j].subtract(_mu));
        } else {
            C[j] = Cj.scalarMultiply(1.d - _r)
                     .add(rxResidual0.outerProduct(x[j].subtract(_mu)));
        }
    }

    // solve A in the following Yule-Walker equation
    // C_j = ∑_{i=1}^{k} A_i C_{j-i} where j = 1..k, C_{-i} = C_i'
    /*
     * /C_1\     /A_1\  /C_0     |C_1'    |C_2'    | .  .  .   |C_{k-1}' \
     * |---|     |---|  |--------+--------+--------+           +---------|
     * |C_2|     |A_2|  |C_1     |C_0     |C_1'    |               .     |
     * |---|     |---|  |--------+--------+--------+               .     |
     * |C_3|  =  |A_3|  |C_2     |C_1     |C_0     |               .     |
     * | . |     | . |  |--------+--------+--------+                     |
     * | . |     | . |  |   .                            .               |
     * | . |     | . |  |   .                            .               |
     * |---|     |---|  |--------+                              +--------|
     * \C_k/     \A_k/  \C_{k-1} | .  .  .                      |C_0     /
     */
    RealMatrix[][] rhs = MatrixUtils.toeplitz(C, k);
    RealMatrix[] lhs = Arrays.copyOfRange(C, 1, k + 1);
    RealMatrix R = MatrixUtils.combinedMatrices(rhs, dims);
    RealMatrix L = MatrixUtils.combinedMatrices(lhs);
    RealMatrix A = MatrixUtils.solve(L, R, false);

    // estimate x
    // \hat{x} = \hat{µ} + ∑_{i=1}^k A_i (x_{t-i} - \hat{µ})
    RealVector x_hat = _mu.copy();
    for (int i = 0; i < k; i++) {
        int offset = i * dims;
        RealMatrix Ai = A.getSubMatrix(offset, offset + dims - 1, 0, dims - 1);
        x_hat = x_hat.add(Ai.operate(xResidual[i + 1]));
    }

    // update model covariance
    // ∑ := (1-r) ∑ + r (x - \hat{x}) (x - \hat{x})'
    RealVector xEstimateResidual = x_t.subtract(x_hat);
    this._sigma = _sigma.scalarMultiply(1.d - _r)
                        .add(xEstimateResidual.mapMultiply(_r).outerProduct(xEstimateResidual));

    return x_hat;
}
 

private void computeADFStatistics() {
	double[] y = diff(ts);
	RealMatrix designMatrix = null;
	int k = lag+1;
	int n = ts.length - 1;
	
	RealMatrix z = MatrixUtils.createRealMatrix(laggedMatrix(y, k)); //has rows length(ts) - 1 - k + 1
	RealVector zcol1 = z.getColumnVector(0); //has length length(ts) - 1 - k + 1
	double[] xt1 = subsetArray(ts, k-1, n-1);  //ts[k:(length(ts) - 1)], has length length(ts) - 1 - k + 1
	double[] trend = sequence(k,n); //trend k:n, has length length(ts) - 1 - k + 1
	if (k > 1) {
		RealMatrix yt1 = z.getSubMatrix(0, ts.length - 1 - k, 1, k-1); //same as z but skips first column
		//build design matrix as cbind(xt1, 1, trend, yt1)
		designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3 + k - 1);
		designMatrix.setColumn(0, xt1);
		designMatrix.setColumn(1, ones(ts.length - 1 - k + 1));
		designMatrix.setColumn(2, trend);
		designMatrix.setSubMatrix(yt1.getData(), 0, 3);
		
	} else {
		//build design matrix as cbind(xt1, 1, tt)
		designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3);
		designMatrix.setColumn(0, xt1);
		designMatrix.setColumn(1, ones(ts.length - 1 - k + 1));
		designMatrix.setColumn(2, trend);
	}
	/*OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
	regression.setNoIntercept(true);
	regression.newSampleData(zcol1.toArray(), designMatrix.getData());
	double[] beta = regression.estimateRegressionParameters();
	double[] sd = regression.estimateRegressionParametersStandardErrors();
	*/
	RidgeRegression regression = new RidgeRegression(designMatrix.getData(), zcol1.toArray());
	regression.updateCoefficients(.0001);
	double[] beta = regression.getCoefficients();
	double[] sd = regression.getStandarderrors();
	
	double t = beta[0] / sd[0];
	if (t <= PVALUE_THRESHOLD) {
		this.needsDiff = true;
	} else {
		this.needsDiff = false;
	}	
}