/** * @param m Input matrix 1. * @param n Input matrix 2. * @return the matrix where the elements of m and n are element-wise multiplied. */ private static RealMatrix times(final RealMatrix m, final RealMatrix n) { final 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] = m.getEntry(r, c) * n.getEntry(r, c); } } return new Array2DRowRealMatrix(d, false); }
final RealVector residuals) { final int nR = jacobian.getRowDimension(); final int nC = jacobian.getColumnDimension(); final RealVector jTr = new ArrayRealVector(nC); jTr.setEntry(j, jTr.getEntry(j) + residuals.getEntry(i) * jacobian.getEntry(i, j)); normal.setEntry(k, l, normal.getEntry(k, l) + jacobian.getEntry(i, k) * jacobian.getEntry(i, l)); normal.setEntry(i, j, normal.getEntry(j, i));
/** * 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; }
/** * Copies a column from m1 to m2. * * @param m1 Source matrix. * @param col1 Source column. * @param m2 Target matrix. * @param col2 Target column. */ private static void copyColumn(final RealMatrix m1, int col1, RealMatrix m2, int col2) { for (int i = 0; i < m1.getRowDimension(); i++) { m2.setEntry(i, col2, m1.getEntry(i, col1)); } }
/** * @param m Input matrix. * @param cols Columns to select. * @return Matrix representing the selected columns. */ private static RealMatrix selectColumns(final RealMatrix m, final int[] cols) { final double[][] d = new double[m.getRowDimension()][cols.length]; for (int r = 0; r < m.getRowDimension(); r++) { for (int c = 0; c < cols.length; c++) { d[r][c] = m.getEntry(r, cols[c]); } } return new Array2DRowRealMatrix(d, false); }
double mean = 0.0; for(int row=0;row<n;row++) { mean += X.getEntry(row, columnId); X.addToEntry(row, columnId, -mean); meanValues.setEntry(columnId, mean); RealMatrix covarianceDD = (X.transpose().multiply(X)).scalarMultiply(1.0/(n-1.0)); RealVector eigenValues = new ArrayRealVector(decomposition.getRealEigenvalues(), false); sqrtEigenValues.setEntry(i, i, FastMath.sqrt(eigenValues.getEntry(i))); components = components.multiply(sqrtEigenValues); double totalVariance = 0.0; for(int i=0;i<d;i++) { totalVariance += eigenValues.getEntry(i); components = components.getSubMatrix(0, components.getRowDimension()-1, 0, maxDimensions-1);
final int n = m.getRowDimension(); if (m.getColumnDimension() != n) { throw new NonSquareMatrixException(m.getRowDimension(), m.getColumnDimension()); final RealMatrix a = m.getSubMatrix(0, splitIndex, 0, splitIndex); final RealMatrix b = m.getSubMatrix(0, splitIndex, splitIndex1, n - 1); final RealMatrix c = m.getSubMatrix(splitIndex1, n - 1, 0, splitIndex); final RealMatrix d = m.getSubMatrix(splitIndex1, n - 1, splitIndex1, n - 1); final RealMatrix aInv = aSolver.getInverse(); final RealMatrix tmp1 = a.subtract(b.multiply(dInv).multiply(c)); final SingularValueDecomposition tmp1Dec = new SingularValueDecomposition(tmp1); final DecompositionSolver tmp1Solver = tmp1Dec.getSolver(); final RealMatrix tmp2 = d.subtract(c.multiply(aInv).multiply(b)); final SingularValueDecomposition tmp2Dec = new SingularValueDecomposition(tmp2); final DecompositionSolver tmp2Solver = tmp2Dec.getSolver(); final RealMatrix result01 = aInv.multiply(b).multiply(result11).scalarMultiply(-1); final RealMatrix result10 = dInv.multiply(c).multiply(result00).scalarMultiply(-1); final RealMatrix result = new Array2DRowRealMatrix(n, n); result.setSubMatrix(result00.getData(), 0, 0); result.setSubMatrix(result01.getData(), 0, splitIndex1); result.setSubMatrix(result10.getData(), splitIndex1, 0); result.setSubMatrix(result11.getData(), splitIndex1, splitIndex1);
Array2DRowRealMatrix H = new Array2DRowRealMatrix(new double[][]{ {(double) dxx, (double) dxy, (double) dxi}, {(double) dxy, (double) dyy, (double) dyi}, RealMatrix H_inv; try { H_inv = new LUDecomposition(H).getSolver().getInverse(); } catch (RuntimeException e) { continue X; double[][] h_inv = H_inv.getData();
RealMatrix Xt = X.transpose(); LUDecomposition lud = new LUDecomposition(Xt.multiply(X)); RealMatrix XtXinv = lud.getSolver().getInverse(); RealVector coefficients = XtXinv.multiply(Xt).operate(Y); thitas.put(Dataframe.COLUMN_NAME_CONSTANT, coefficients.getEntry(0)); for(Map.Entry<Object, Integer> entry : featureIds.entrySet()) { Object feature = entry.getKey(); Integer featureId = entry.getValue(); thitas.put(feature, coefficients.getEntry(featureId)); for(double v : X.operate(coefficients).subtract(Y).toArray()) { SSE += v*v; RealMatrix SE = XtXinv.scalarMultiply(MSE); double error = SE.getEntry(i, i); Object feature = idsFeatures.get(i); if(error<=0.0) {
/** {@inheritDoc} */ public RealVector getSigma(double covarianceSingularityThreshold) { final RealMatrix cov = this.getCovariances(covarianceSingularityThreshold); final int nC = cov.getColumnDimension(); final RealVector sig = new ArrayRealVector(nC); for (int i = 0; i < nC; ++i) { sig.setEntry(i, FastMath.sqrt(cov.getEntry(i,i))); } return sig; }
residualsWeights[i] = weightMatrix.getEntry(i, i); for (int i = 0; i < nR; ++i) { final double[] grad = weightedJacobian.getRow(i); final double weight = residualsWeights[i]; final double residual = currentResiduals[i]; new LUDecomposition(mA).getSolver() : new QRDecomposition(mA).getSolver(); final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray();
double negccov = 0; if (ccov1 + ccovmu > 0) { final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu)) .scalarMultiply(1 / sigma); // mu difference vectors final RealMatrix roneu = pc.multiply(pc.transpose()) .scalarMultiply(ccov1); // rank one update RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu)); RealMatrix arnorms = sqrt(sumRows(square(arzneg))); final int[] idxnorms = sortedIndices(arnorms.getRow(0)); final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms); final int[] idxReverse = reverse(idxnorms); square(arnormsInv).multiply(weights).getEntry(0, 0); if (negccov > negcovMax) { negccov = negcovMax; final RealMatrix artmp = BD.multiply(arzneg); final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose()); oldFac += negalphaold * negccov; C = C.scalarMultiply(oldFac) .add(roneu) // regard old matrix .add(arpos.scalarMultiply( // plus rank one update ccovmu + (1 - negalphaold) * negccov) // plus rank mu update .multiply(times(repmat(weights, 1, dimension), arpos.transpose()))) .subtract(Cneg.scalarMultiply(negccov)); } else {
transitionMatrixT = transitionMatrix.transpose(); measurementMatrixT = measurementMatrix.transpose(); stateEstimation = new ArrayRealVector(transitionMatrix.getColumnDimension()); } else { stateEstimation = processModel.getInitialStateEstimate(); if (transitionMatrix.getColumnDimension() != stateEstimation.getDimension()) { throw new DimensionMismatchException(transitionMatrix.getColumnDimension(), stateEstimation.getDimension()); errorCovariance = processNoise.copy(); } else { errorCovariance = processModel.getInitialErrorCovariance(); if (!transitionMatrix.isSquare()) { throw new NonSquareMatrixException( transitionMatrix.getRowDimension(), transitionMatrix.getColumnDimension()); controlMatrix.getRowDimension() > 0 && controlMatrix.getColumnDimension() > 0 && controlMatrix.getRowDimension() != transitionMatrix.getRowDimension()) { throw new MatrixDimensionMismatchException(controlMatrix.getRowDimension(), controlMatrix.getColumnDimension(), transitionMatrix.getRowDimension(), controlMatrix.getColumnDimension());
/** * Computes the square-root of the weight matrix. * * @param m Symmetric, positive-definite (weight) matrix. * @return the square-root of the weight matrix. */ private static RealMatrix squareRoot(final RealMatrix m) { if (m instanceof DiagonalMatrix) { final int dim = m.getRowDimension(); final RealMatrix sqrtM = new DiagonalMatrix(dim); for (int i = 0; i < dim; i++) { sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i))); } return sqrtM; } else { final EigenDecomposition dec = new EigenDecomposition(m); return dec.getSquareRoot(); } }
if (z.getDimension() != measurementMatrix.getRowDimension()) { throw new DimensionMismatchException(z.getDimension(), measurementMatrix.getRowDimension()); RealMatrix s = measurementMatrix.multiply(errorCovariance) .multiply(measurementMatrixT) .add(measurementModel.getMeasurementNoise()); RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation)); .solve(measurementMatrix.multiply(errorCovariance.transpose())) .transpose(); stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance);
public static double sumWalks(RealMatrix m, RealVector s, RealVector t) { int n = m.getColumnDimension(); if (n != m.getRowDimension()) { throw new InputMismatchException("sum walks can only be computed on square matrices"); } LUDecomposition d = new LUDecomposition(MatrixUtils.createRealIdentityMatrix(n).subtract(m)); RealMatrix sum = d.getSolver().getInverse(); if (s == null) s = new ArrayRealVector(n, 1.0); if (t == null) t = new ArrayRealVector(n, 1.0); return sum.preMultiply(s).dotProduct(t); }
/** * <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()); }
/** * Calculates the QR-decomposition of the given matrix. * * @param matrix The matrix to decompose. * @param threshold Singularity threshold. */ public QRDecomposition(RealMatrix matrix, double threshold) { this.threshold = threshold; final int m = matrix.getRowDimension(); final int n = matrix.getColumnDimension(); qrt = matrix.transpose().getData(); rDiag = new double[FastMath.min(m, n)]; cachedQ = null; cachedQT = null; cachedR = null; cachedH = null; decompose(qrt); }
/** * 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()); }
if (b.getRowDimension() != m) { throw new DimensionMismatchException(b.getRowDimension(), m); final int nColB = b.getColumnDimension(); final double[][] bp = new double[m][nColB]; final double[] tmpCol = new double[m]; for (int k = 0; k < nColB; ++k) { for (int i = 0; i < m; ++i) { tmpCol[i] = b.getEntry(i, k); bp[i][k] = 0; final double[] vData = v.getDataRef(); double s = 0; for (int j = 0; j < m; ++j) { s += v.getEntry(j) * tmpCol[j]; return new Array2DRowRealMatrix(bp, false);