List of usage examples for org.apache.commons.math3.linear QRDecomposition getSolver
public DecompositionSolver getSolver()
From source file:com.ibm.bi.dml.runtime.matrix.data.LibCommonsMath.java
/** * Function to compute matrix inverse via matrix decomposition. * /*from w ww.j a v a 2 s .c om*/ * @param in * @return * @throws DMLRuntimeException */ private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix in) throws DMLRuntimeException { if (!in.isSquare()) throw new DMLRuntimeException("Input to inv() must be square matrix -- given: a " + in.getRowDimension() + "x" + in.getColumnDimension() + " matrix."); QRDecomposition qrdecompose = new QRDecomposition(in); DecompositionSolver solver = qrdecompose.getSolver(); RealMatrix inverseMatrix = solver.getInverse(); MatrixBlock inverse = DataConverter.convertToMatrixBlock(inverseMatrix.getData()); return inverse; }
From source file:com.ibm.bi.dml.runtime.matrix.data.LibCommonsMath.java
/** * Function to solve a given system of equations. * /* ww w. ja v a 2s.c o m*/ * @param in1 * @param in2 * @return * @throws DMLRuntimeException */ private static MatrixBlock computeSolve(MatrixObject in1, MatrixObject in2) throws DMLRuntimeException { Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in1); Array2DRowRealMatrix vectorInput = DataConverter.convertToArray2DRowRealMatrix(in2); /*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput); DecompositionSolver lusolver = ludecompose.getSolver(); RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/ // Setup a solver based on QR Decomposition QRDecomposition qrdecompose = new QRDecomposition(matrixInput); DecompositionSolver solver = qrdecompose.getSolver(); // Invoke solve RealMatrix solutionMatrix = solver.solve(vectorInput); MatrixBlock solution = DataConverter.convertToMatrixBlock(solutionMatrix.getData()); return solution; }
From source file:com.joptimizer.util.Utils.java
/** * @see http://en.wikipedia.org/wiki/Machine_epsilon#Approximation_using_Java *//* ww w.jav a2 s. c o m*/ // public static final float calculateFloatMachineEpsilon() { // float eps = 1.0f; // do { // eps /= 2.0f; // } while ((float) (1.0 + (eps / 2.0)) != 1.0); // Log.d(MainActivity.JOPTIMIZER_LOGTAG,"Calculated float machine epsilon: " + eps); // return eps; // } public static RealMatrix squareMatrixInverse(RealMatrix M) throws SingularMatrixException { if (!M.isSquare()) { throw new IllegalArgumentException("Not square matrix!"); } // try commons-math cholesky try { CholeskyDecomposition cd = new CholeskyDecomposition(M); return cd.getSolver().getInverse(); } catch (SingularMatrixException e) { throw e; } catch (Exception e) { Log.d(MainActivity.JOPTIMIZER_LOGTAG, e.getMessage()); } // try joptimizer cholesky try { CholeskyFactorization cd = new CholeskyFactorization(M.getData()); double[][] MInv = cd.getInverse(); return MatrixUtils.createRealMatrix(MInv); } catch (Exception e) { Log.d(MainActivity.JOPTIMIZER_LOGTAG, e.getMessage()); } // try LU try { LUDecomposition ld = new LUDecomposition(M); return ld.getSolver().getInverse(); } catch (SingularMatrixException e) { throw e; } catch (Exception e) { Log.d(MainActivity.JOPTIMIZER_LOGTAG, e.getMessage()); } // try QR try { QRDecomposition qr = new org.apache.commons.math3.linear.QRDecomposition(M); return qr.getSolver().getInverse(); } catch (SingularMatrixException e) { throw e; } catch (Exception e) { Log.d(MainActivity.JOPTIMIZER_LOGTAG, e.getMessage()); } return null; }
From source file:com.anhth12.lambda.common.math.LinearSystemSolver.java
public boolean isNonSingular(RealMatrix M) { QRDecomposition decomposition = new RRQRDecomposition(M, SINGULARITY_THRESHOLD); DecompositionSolver solver = decomposition.getSolver(); return solver.isNonSingular(); }
From source file:com.cloudera.oryx.common.math.CommonsMathLinearSystemSolver.java
@Override public boolean isNonSingular(RealMatrix M) { QRDecomposition decomposition = new RRQRDecomposition(M, SINGULARITY_THRESHOLD); DecompositionSolver solver = decomposition.getSolver(); return solver.isNonSingular(); }
From source file:com.cloudera.oryx.common.math.LinearSystemSolver.java
public boolean isNonSingular(RealMatrix M) { QRDecomposition decomposition = new RRQRDecomposition(M, SINGULARITY_THRESHOLD); DecompositionSolver solver = decomposition.getSolver(); return solver.isNonSingular(); }
From source file:com.opengamma.strata.math.impl.linearalgebra.QRDecompositionCommonsResult.java
/** * @param qr The result of the QR decomposition, not null *///from w ww . j av a2 s.co m public QRDecompositionCommonsResult(QRDecomposition qr) { ArgChecker.notNull(qr, "qr"); _q = CommonsMathWrapper.unwrap(qr.getQ()); _r = CommonsMathWrapper.unwrap(qr.getR()); _qTranspose = _q.transpose(); _solver = qr.getSolver(); }
From source file:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Calculates the leverages of data points for least squares fitting (assuming equal variances). * * @param indVarValues The values of the independent variable used for the fitting. * @return a RealVector containing a leverage value for each independent variable value. *///from w ww .j a va 2s . co m protected RealVector calculateLeverages(RealVector indVarValues) { RealMatrix indVarMatrix = null; if (this.noIntercept) { indVarMatrix = new Array2DRowRealMatrix(indVarValues.getDimension(), 1); } else { indVarMatrix = new Array2DRowRealMatrix(indVarValues.getDimension(), 2); } indVarMatrix.setColumnVector(0, indVarValues); if (!this.noIntercept) { indVarMatrix.setColumnVector(1, indVarValues.mapMultiply(0).mapAdd(1)); } RealVector leverages = new ArrayRealVector(indVarValues.getDimension()); QRDecomposition xQR = new QRDecomposition(indVarMatrix); RealMatrix xR = xQR.getR(); int smallerDim = xR.getRowDimension() < xR.getColumnDimension() ? xR.getRowDimension() : xR.getColumnDimension(); RealMatrix xRSq = xR.getSubMatrix(0, smallerDim - 1, 0, smallerDim - 1); QRDecomposition xRQR = new QRDecomposition(xRSq); RealMatrix xRInv = xRQR.getSolver().getInverse(); RealMatrix xxRInv = indVarMatrix.multiply(xRInv); for (int i = 0; i < indVarValues.getDimension(); i++) { double sum = 0; for (int j = 0; j < xxRInv.getColumnDimension(); j++) { sum += Math.pow(xxRInv.getEntry(i, j), 2); } leverages.setEntry(i, sum); } return leverages; }
From source file:com.google.location.lbs.gnss.gps.pseudorange.UserPositionVelocityWeightedLeastSquare.java
/** * Least square solution to calculate the user position given the navigation message, pseudorange * and accumulated delta range measurements. Also calculates user velocity non-iteratively from * Least square position solution.//from w w w . j a va2 s . c o m * * <p>The method fills the user position and velocity in ECEF coordinates and receiver clock * offset in meters and clock offset rate in meters per second. * * <p>One can choose between no smoothing, using the carrier phase measurements (accumulated delta * range) or the doppler measurements (pseudorange rate) for smoothing the pseudorange. The * smoothing is applied only if time has changed below a specific threshold since last invocation. * * <p>Source for least squares: * * <ul> * <li>http://www.u-blox.com/images/downloads/Product_Docs/GPS_Compendium%28GPS-X-02007%29.pdf * page 81 - 85 * <li>Parkinson, B.W., Spilker Jr., J.J.: Global positioning system: theory and applications * page 412 - 414 * </ul> * * <p>Sources for smoothing pseudorange with carrier phase measurements: * * <ul> * <li>Satellite Communications and Navigation Systems book, page 424, * <li>Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems, page 388, * 389. * </ul> * * <p>The function does not modify the smoothed measurement list {@code * immutableSmoothedSatellitesToReceiverMeasurements} * * @param navMessageProto parameters of the navigation message * @param usefulSatellitesToReceiverMeasurements Map of useful satellite PRN to {@link * GpsMeasurementWithRangeAndUncertainty} containing receiver measurements for computing the * position solution. * @param receiverGPSTowAtReceptionSeconds Receiver estimate of GPS time of week (seconds) * @param receiverGPSWeek Receiver estimate of GPS week (0-1024+) * @param dayOfYear1To366 The day of the year between 1 and 366 * @param positionVelocitySolutionECEF Solution array of the following format: * [0-2] xyz solution of user. * [3] clock bias of user. * [4-6] velocity of user. * [7] clock bias rate of user. * @param positionVelocityUncertaintyEnu Uncertainty of calculated position and velocity solution * in meters and mps local ENU system. Array has the following format: * [0-2] Enu uncertainty of position solution in meters * [3-5] Enu uncertainty of velocity solution in meters per second. * @param pseudorangeResidualMeters The pseudorange residual corrected by subtracting expected * psudorange calculated with the use clock bias of the highest elevation satellites. */ public void calculateUserPositionVelocityLeastSquare(GpsNavMessageProto navMessageProto, List<GpsMeasurementWithRangeAndUncertainty> usefulSatellitesToReceiverMeasurements, double receiverGPSTowAtReceptionSeconds, int receiverGPSWeek, int dayOfYear1To366, double[] positionVelocitySolutionECEF, double[] positionVelocityUncertaintyEnu, double[] pseudorangeResidualMeters) throws Exception { // Use PseudorangeSmoother to smooth the pseudorange according to: Satellite Communications and // Navigation Systems book, page 424 and Principles of GNSS, Inertial, and Multisensor // Integrated Navigation Systems, page 388, 389. double[] deltaPositionMeters; List<GpsMeasurementWithRangeAndUncertainty> immutableSmoothedSatellitesToReceiverMeasurements = pseudorangeSmoother .updatePseudorangeSmoothingResult( Collections.unmodifiableList(usefulSatellitesToReceiverMeasurements)); List<GpsMeasurementWithRangeAndUncertainty> mutableSmoothedSatellitesToReceiverMeasurements = Lists .newArrayList(immutableSmoothedSatellitesToReceiverMeasurements); int numberOfUsefulSatellites = getNumberOfUsefulSatellites(mutableSmoothedSatellitesToReceiverMeasurements); // Least square position solution is supported only if 4 or more satellites visible Preconditions.checkArgument(numberOfUsefulSatellites >= MINIMUM_NUMER_OF_SATELLITES, "At least 4 satellites have to be visible... Only 3D mode is supported..."); boolean repeatLeastSquare = false; SatellitesPositionPseudorangesResidualAndCovarianceMatrix satPosPseudorangeResidualAndWeight; boolean isFirstWLS = true; do { // Calculate satellites' positions, measurement residuals per visible satellite and // weight matrix for the iterative least square boolean doAtmosphericCorrections = false; satPosPseudorangeResidualAndWeight = calculateSatPosAndPseudorangeResidual(navMessageProto, mutableSmoothedSatellitesToReceiverMeasurements, receiverGPSTowAtReceptionSeconds, receiverGPSWeek, dayOfYear1To366, positionVelocitySolutionECEF, doAtmosphericCorrections); // Calculate the geometry matrix according to "Global Positioning System: Theory and // Applications", Parkinson and Spilker page 413 RealMatrix covarianceMatrixM2 = new Array2DRowRealMatrix( satPosPseudorangeResidualAndWeight.covarianceMatrixMetersSquare); geometryMatrix = new Array2DRowRealMatrix(calculateGeometryMatrix( satPosPseudorangeResidualAndWeight.satellitesPositionsMeters, positionVelocitySolutionECEF)); RealMatrix weightedGeometryMatrix; RealMatrix weightMatrixMetersMinus2 = null; // Apply weighted least square only if the covariance matrix is not singular (has a non-zero // determinant), otherwise apply ordinary least square. The reason is to ignore reported // signal to noise ratios by the receiver that can lead to such singularities LUDecomposition ludCovMatrixM2 = new LUDecomposition(covarianceMatrixM2); double det = ludCovMatrixM2.getDeterminant(); if (det <= DOUBLE_ROUND_OFF_TOLERANCE) { // Do not weight the geometry matrix if covariance matrix is singular. weightedGeometryMatrix = geometryMatrix; } else { weightMatrixMetersMinus2 = ludCovMatrixM2.getSolver().getInverse(); RealMatrix hMatrix = calculateHMatrix(weightMatrixMetersMinus2, geometryMatrix); weightedGeometryMatrix = hMatrix.multiply(geometryMatrix.transpose()) .multiply(weightMatrixMetersMinus2); } // Equation 9 page 413 from "Global Positioning System: Theory and Applicaitons", Parkinson // and Spilker deltaPositionMeters = GpsMathOperations.matrixByColVectMultiplication(weightedGeometryMatrix.getData(), satPosPseudorangeResidualAndWeight.pseudorangeResidualsMeters); // Apply corrections to the position estimate positionVelocitySolutionECEF[0] += deltaPositionMeters[0]; positionVelocitySolutionECEF[1] += deltaPositionMeters[1]; positionVelocitySolutionECEF[2] += deltaPositionMeters[2]; positionVelocitySolutionECEF[3] += deltaPositionMeters[3]; // Iterate applying corrections to the position solution until correction is below threshold satPosPseudorangeResidualAndWeight = applyWeightedLeastSquare(navMessageProto, mutableSmoothedSatellitesToReceiverMeasurements, receiverGPSTowAtReceptionSeconds, receiverGPSWeek, dayOfYear1To366, positionVelocitySolutionECEF, deltaPositionMeters, doAtmosphericCorrections, satPosPseudorangeResidualAndWeight, weightMatrixMetersMinus2); // We use the first WLS iteration results and correct them based on the ground truth position // and using a clock error computed from high elevation satellites. The first iteration is // used before satellite with high residuals being removed. if (isFirstWLS && truthLocationForCorrectedResidualComputationEcef != null) { // Snapshot the information needed before high residual satellites are removed System.arraycopy( ResidualCorrectionCalculator.calculateCorrectedResiduals(satPosPseudorangeResidualAndWeight, positionVelocitySolutionECEF.clone(), truthLocationForCorrectedResidualComputationEcef), 0 /*source starting pos*/, pseudorangeResidualMeters, 0 /*destination starting pos*/, GpsNavigationMessageStore.MAX_NUMBER_OF_SATELLITES /*length of elements*/); isFirstWLS = false; } repeatLeastSquare = false; int satsWithResidualBelowThreshold = satPosPseudorangeResidualAndWeight.pseudorangeResidualsMeters.length; // remove satellites that have residuals above RESIDUAL_TO_REPEAT_LEAST_SQUARE_METERS as they // worsen the position solution accuracy. If any satellite is removed, repeat the least square repeatLeastSquare = removeHighResidualSats(mutableSmoothedSatellitesToReceiverMeasurements, repeatLeastSquare, satPosPseudorangeResidualAndWeight, satsWithResidualBelowThreshold); } while (repeatLeastSquare); calculateGeoidMeters = false; // The computed ECEF position will be used next to compute the user velocity. // we calculate and fill in the user velocity solutions based on following equation: // Weight Matrix * GeometryMatrix * User Velocity Vector // = Weight Matrix * deltaPseudoRangeRateWeightedMps // Reference: Pratap Misra and Per Enge // "Global Positioning System: Signals, Measurements, and Performance" Page 218. // Get the number of satellite used in Geometry Matrix numberOfUsefulSatellites = geometryMatrix.getRowDimension(); RealMatrix rangeRateMps = new Array2DRowRealMatrix(numberOfUsefulSatellites, 1); RealMatrix deltaPseudoRangeRateMps = new Array2DRowRealMatrix(numberOfUsefulSatellites, 1); RealMatrix pseudorangeRateWeight = new Array2DRowRealMatrix(numberOfUsefulSatellites, numberOfUsefulSatellites); // Correct the receiver time of week with the estimated receiver clock bias receiverGPSTowAtReceptionSeconds = receiverGPSTowAtReceptionSeconds - positionVelocitySolutionECEF[3] / SPEED_OF_LIGHT_MPS; int measurementCount = 0; // Calculate range rates for (int i = 0; i < GpsNavigationMessageStore.MAX_NUMBER_OF_SATELLITES; i++) { if (mutableSmoothedSatellitesToReceiverMeasurements.get(i) != null) { GpsEphemerisProto ephemeridesProto = getEphemerisForSatellite(navMessageProto, i + 1); double pseudorangeMeasurementMeters = mutableSmoothedSatellitesToReceiverMeasurements .get(i).pseudorangeMeters; GpsTimeOfWeekAndWeekNumber correctedTowAndWeek = calculateCorrectedTransmitTowAndWeek( ephemeridesProto, receiverGPSTowAtReceptionSeconds, receiverGPSWeek, pseudorangeMeasurementMeters); // Calculate satellite velocity PositionAndVelocity satPosECEFMetersVelocityMPS = SatellitePositionCalculator .calculateSatellitePositionAndVelocityFromEphemeris(ephemeridesProto, correctedTowAndWeek.gpsTimeOfWeekSeconds, correctedTowAndWeek.weekNumber, positionVelocitySolutionECEF[0], positionVelocitySolutionECEF[1], positionVelocitySolutionECEF[2]); // Calculate satellite clock error rate double satelliteClockErrorRateMps = SatelliteClockCorrectionCalculator .calculateSatClockCorrErrorRate(ephemeridesProto, correctedTowAndWeek.gpsTimeOfWeekSeconds, correctedTowAndWeek.weekNumber); // Fill in range rates. range rate = satellite velocity (dot product) line-of-sight vector rangeRateMps.setEntry(measurementCount, 0, -1 * (satPosECEFMetersVelocityMPS.velocityXMetersPerSec * geometryMatrix.getEntry(measurementCount, 0) + satPosECEFMetersVelocityMPS.velocityYMetersPerSec * geometryMatrix.getEntry(measurementCount, 1) + satPosECEFMetersVelocityMPS.velocityZMetersPerSec * geometryMatrix.getEntry(measurementCount, 2))); deltaPseudoRangeRateMps.setEntry(measurementCount, 0, mutableSmoothedSatellitesToReceiverMeasurements.get(i).pseudorangeRateMps - rangeRateMps.getEntry(measurementCount, 0) + satelliteClockErrorRateMps - positionVelocitySolutionECEF[7]); // Calculate the velocity weight matrix by using 1 / square(Pseudorangerate Uncertainty) // along the diagonal pseudorangeRateWeight.setEntry(measurementCount, measurementCount, 1 / (mutableSmoothedSatellitesToReceiverMeasurements.get(i).pseudorangeRateUncertaintyMps * mutableSmoothedSatellitesToReceiverMeasurements .get(i).pseudorangeRateUncertaintyMps)); measurementCount++; } } RealMatrix weightedGeoMatrix = pseudorangeRateWeight.multiply(geometryMatrix); RealMatrix deltaPseudoRangeRateWeightedMps = pseudorangeRateWeight.multiply(deltaPseudoRangeRateMps); QRDecomposition qrdWeightedGeoMatrix = new QRDecomposition(weightedGeoMatrix); RealMatrix velocityMps = qrdWeightedGeoMatrix.getSolver().solve(deltaPseudoRangeRateWeightedMps); positionVelocitySolutionECEF[4] = velocityMps.getEntry(0, 0); positionVelocitySolutionECEF[5] = velocityMps.getEntry(1, 0); positionVelocitySolutionECEF[6] = velocityMps.getEntry(2, 0); positionVelocitySolutionECEF[7] = velocityMps.getEntry(3, 0); RealMatrix pseudorangeWeight = new LUDecomposition( new Array2DRowRealMatrix(satPosPseudorangeResidualAndWeight.covarianceMatrixMetersSquare)) .getSolver().getInverse(); // Calculate and store the uncertainties of position and velocity in local ENU system in meters // and meters per second. double[] pvUncertainty = calculatePositionVelocityUncertaintyEnu(pseudorangeRateWeight, pseudorangeWeight, positionVelocitySolutionECEF); System.arraycopy(pvUncertainty, 0 /*source starting pos*/, positionVelocityUncertaintyEnu, 0 /*destination starting pos*/, 6 /*length of elements*/); }
From source file:org.apache.sysml.runtime.matrix.data.LibCommonsMath.java
/** * Function to solve a given system of equations. * // w ww . java2 s .c o m * @param in1 matrix object 1 * @param in2 matrix object 2 * @return matrix block * @throws DMLRuntimeException if DMLRuntimeException occurs */ private static MatrixBlock computeSolve(MatrixObject in1, MatrixObject in2) throws DMLRuntimeException { Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in1); Array2DRowRealMatrix vectorInput = DataConverter.convertToArray2DRowRealMatrix(in2); /*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput); DecompositionSolver lusolver = ludecompose.getSolver(); RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/ // Setup a solver based on QR Decomposition QRDecomposition qrdecompose = new QRDecomposition(matrixInput); DecompositionSolver solver = qrdecompose.getSolver(); // Invoke solve RealMatrix solutionMatrix = solver.solve(vectorInput); return DataConverter.convertToMatrixBlock(solutionMatrix.getData()); }