/** {@inheritDoc} */ public RealMatrix preMultiply(final RealMatrix m) throws DimensionMismatchException { return m.multiply(this); }
/** * Solve the linear equation A × X = B in least square sense. * <p> * The m×n matrix A may not be square, the solution X is such that * ||A × X - B|| is minimal. * </p> * * @param b Right-hand side of the equation A × X = B * @return a matrix X that minimizes the two norm of A × X - B * @throws org.apache.commons.math3.exception.DimensionMismatchException * if the matrices dimensions do not match. */ public RealMatrix solve(final RealMatrix b) { return pseudoInverse.multiply(b); }
/** {@inheritDoc} */ public RealMatrix getJacobian() { return weightSqrt.multiply(this.unweighted.getJacobian()); }
/** {@inheritDoc} */ public RealMatrix solve(RealMatrix b) { return p.multiply(upper.solve(b)); }
/** * Computes the weighted Jacobian matrix. * * @param params Model parameters at which to compute the Jacobian. * @return the weighted Jacobian: W<sup>1/2</sup> J. * @throws DimensionMismatchException if the Jacobian dimension does not * match problem dimension. */ protected RealMatrix computeWeightedJacobian(double[] params) { return weightMatrixSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params))); }
/** * 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()); }
/** * Build a solver from decomposed matrix. * * @param singularValues Singular values. * @param uT U<sup>T</sup> matrix of the decomposition. * @param v V matrix of the decomposition. * @param nonSingular Singularity indicator. * @param tol tolerance for singular values */ private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v, final boolean nonSingular, final double tol) { final double[][] suT = uT.getData(); for (int i = 0; i < singularValues.length; ++i) { final double a; if (singularValues[i] > tol) { a = 1 / singularValues[i]; } else { a = 0; } final double[] suTi = suT[i]; for (int j = 0; j < suTi.length; ++j) { suTi[j] *= a; } } pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false)); this.nonSingular = nonSingular; }
/** * Calculates the variance on the beta. * <pre> * Var(b)=(X' Omega^-1 X)^-1 * </pre> * @return The beta variance matrix */ @Override protected RealMatrix calculateBetaVariance() { RealMatrix OI = getOmegaInverse(); RealMatrix XTOIX = getX().transpose().multiply(OI).multiply(getX()); return new LUDecomposition(XTOIX).getSolver().getInverse(); }
/** * Computes the square-root of the matrix. * This implementation assumes that the matrix is symmetric and positive * definite. * * @return the square-root of the matrix. * @throws MathUnsupportedOperationException if the matrix is not * symmetric or not positive definite. * @since 3.1 */ public RealMatrix getSquareRoot() { if (!isSymmetric) { throw new MathUnsupportedOperationException(); } final double[] sqrtEigenValues = new double[realEigenvalues.length]; for (int i = 0; i < realEigenvalues.length; i++) { final double eigen = realEigenvalues[i]; if (eigen <= 0) { throw new MathUnsupportedOperationException(); } sqrtEigenValues[i] = FastMath.sqrt(eigen); } final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues); final RealMatrix v = getV(); final RealMatrix vT = getVT(); return v.multiply(sqrtEigen).multiply(vT); }
/** {@inheritDoc} */ public RealMatrix getCovariances(double threshold) { // Set up the Jacobian. final RealMatrix j = this.getJacobian(); // Compute transpose(J)J. final RealMatrix jTj = j.transpose().multiply(j); // Compute the covariances matrix. final DecompositionSolver solver = new QRDecomposition(jTj, threshold).getSolver(); return solver.getInverse(); }
/** * Predict the internal state estimation one time step ahead. * * @param u * the control vector * @throws DimensionMismatchException * if the dimension of the control vector does not match */ public void predict(final RealVector u) throws DimensionMismatchException { // sanity checks if (u != null && u.getDimension() != controlMatrix.getColumnDimension()) { throw new DimensionMismatchException(u.getDimension(), controlMatrix.getColumnDimension()); } // project the state estimation ahead (a priori state) // xHat(k)- = A * xHat(k-1) + B * u(k-1) stateEstimation = transitionMatrix.operate(stateEstimation); // add control input if it is available if (u != null) { stateEstimation = stateEstimation.add(controlMatrix.operate(u)); } // project the error covariance ahead // P(k)- = A * P(k-1) * A' + Q errorCovariance = transitionMatrix.multiply(errorCovariance) .multiply(transitionMatrixT) .add(processModel.getProcessNoise()); }
/** * Get the covariance matrix of the optimized parameters. * <br/> * Note that this operation involves the inversion of the * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the * Jacobian matrix. * The {@code threshold} parameter is a way for the caller to specify * that the result of this computation should be considered meaningless, * and thus trigger an exception. * * @param params Model parameters. * @param threshold Singularity threshold. * @return the covariance matrix. * @throws org.apache.commons.math3.linear.SingularMatrixException * if the covariance matrix cannot be computed (singular problem). * @since 3.1 */ public double[][] computeCovariances(double[] params, double threshold) { // Set up the Jacobian. final RealMatrix j = computeWeightedJacobian(params); // Compute transpose(J)J. final RealMatrix jTj = j.transpose().multiply(j); // Compute the covariances matrix. final DecompositionSolver solver = new QRDecomposition(jTj, threshold).getSolver(); return solver.getInverse().getData(); }
/** * Get the covariance matrix of the optimized parameters. * <br/> * Note that this operation involves the inversion of the * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the * Jacobian matrix. * The {@code threshold} parameter is a way for the caller to specify * that the result of this computation should be considered meaningless, * and thus trigger an exception. * * @param params Model parameters. * @param threshold Singularity threshold. * @return the covariance matrix. * @throws org.apache.commons.math3.linear.SingularMatrixException * if the covariance matrix cannot be computed (singular problem). */ public double[][] computeCovariances(double[] params, double threshold) { // Set up the Jacobian. final RealMatrix j = computeWeightedJacobian(params); // Compute transpose(J)J. final RealMatrix jTj = j.transpose().multiply(j); // Compute the covariances matrix. final DecompositionSolver solver = new QRDecomposition(jTj, threshold).getSolver(); return solver.getInverse().getData(); }
final RealMatrix X = matrixDataset.getX().multiply(components);
/** * <p>Calculates the variance-covariance matrix of the regression parameters. * </p> * <p>Var(b) = (X<sup>T</sup>X)<sup>-1</sup> * </p> * <p>Uses QR decomposition to reduce (X<sup>T</sup>X)<sup>-1</sup> * to (R<sup>T</sup>R)<sup>-1</sup>, with only the top p rows of * R included, where p = the length of the beta vector.</p> * * <p>Data for the model must have been successfully loaded using one of * the {@code newSampleData} methods before invoking this method; otherwise * a {@code NullPointerException} will be thrown.</p> * * @return The beta variance-covariance matrix * @throws org.apache.commons.math3.linear.SingularMatrixException if the design matrix is singular * @throws NullPointerException if the data for the model have not been loaded */ @Override protected RealMatrix calculateBetaVariance() { int p = getX().getColumnDimension(); RealMatrix Raug = qr.getR().getSubMatrix(0, p - 1 , 0, p - 1); RealMatrix Rinv = new LUDecomposition(Raug).getSolver().getInverse(); return Rinv.multiply(Rinv.transpose()); }
RealMatrix rmResult = rmFirst.multiply(rmSecond);
/** * Update of the evolution paths ps and pc. * * @param zmean Weighted row matrix of the gaussian random numbers generating * the current offspring. * @param xold xmean matrix of the previous generation. * @return hsig flag indicating a small correction. */ private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) { ps = ps.scalarMultiply(1 - cs).add( B.multiply(zmean).scalarMultiply(FastMath.sqrt(cs * (2 - cs) * mueff))); normps = ps.getFrobeniusNorm(); final boolean hsig = normps / FastMath.sqrt(1 - FastMath.pow(1 - cs, 2 * iterations)) / chiN < 1.4 + 2 / ((double) dimension + 1); pc = pc.scalarMultiply(1 - cc); if (hsig) { pc = pc.add(xmean.subtract(xold).scalarMultiply(FastMath.sqrt(cc * (2 - cc) * mueff) / sigma)); } return hsig; }
/** * Update of the evolution paths ps and pc. * * @param zmean Weighted row matrix of the gaussian random numbers generating * the current offspring. * @param xold xmean matrix of the previous generation. * @return hsig flag indicating a small correction. */ private boolean updateEvolutionPaths(RealMatrix zmean, RealMatrix xold) { ps = ps.scalarMultiply(1 - cs).add( B.multiply(zmean).scalarMultiply( FastMath.sqrt(cs * (2 - cs) * mueff))); normps = ps.getFrobeniusNorm(); final boolean hsig = normps / FastMath.sqrt(1 - FastMath.pow(1 - cs, 2 * iterations)) / chiN < 1.4 + 2 / ((double) dimension + 1); pc = pc.scalarMultiply(1 - cc); if (hsig) { pc = pc.add(xmean.subtract(xold).scalarMultiply(FastMath.sqrt(cc * (2 - cc) * mueff) / sigma)); } return hsig; }
/** * Update of the covariance matrix C for diagonalOnly > 0 * * @param hsig Flag indicating a small correction. * @param bestArz Fitness-sorted matrix of the gaussian random values of the * current offspring. */ private void updateCovarianceDiagonalOnly(boolean hsig, final RealMatrix bestArz) { // minor correction if hsig==false double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc); oldFac += 1 - ccov1Sep - ccovmuSep; diagC = diagC.scalarMultiply(oldFac) // regard old matrix .add(square(pc).scalarMultiply(ccov1Sep)) // plus rank one update .add((times(diagC, square(bestArz).multiply(weights))) // plus rank mu update .scalarMultiply(ccovmuSep)); diagD = sqrt(diagC); // replaces eig(C) if (diagonalOnly > 1 && iterations > diagonalOnly) { // full covariance matrix from now on diagonalOnly = 0; B = eye(dimension, dimension); BD = diag(diagD); C = diag(diagC); } }
/** * Update of the covariance matrix C for diagonalOnly > 0 * * @param hsig Flag indicating a small correction. * @param bestArz Fitness-sorted matrix of the gaussian random values of the * current offspring. */ private void updateCovarianceDiagonalOnly(boolean hsig, final RealMatrix bestArz) { // minor correction if hsig==false double oldFac = hsig ? 0 : ccov1Sep * cc * (2 - cc); oldFac += 1 - ccov1Sep - ccovmuSep; diagC = diagC.scalarMultiply(oldFac) // regard old matrix .add(square(pc).scalarMultiply(ccov1Sep)) // plus rank one update .add((times(diagC, square(bestArz).multiply(weights))) // plus rank mu update .scalarMultiply(ccovmuSep)); diagD = sqrt(diagC); // replaces eig(C) if (diagonalOnly > 1 && iterations > diagonalOnly) { // full covariance matrix from now on diagonalOnly = 0; B = eye(dimension, dimension); BD = diag(diagD); C = diag(diagC); } }