use of org.apache.commons.math3.linear.SingularMatrixException in project GDSC-SMLM by aherbert.
the class ApacheLVMFitter method computeFit.
public FitStatus computeFit(double[] y, final double[] y_fit, double[] a, double[] a_dev) {
int n = y.length;
try {
// Different convergence thresholds seem to have no effect on the resulting fit, only the number of
// iterations for convergence
final double initialStepBoundFactor = 100;
final double costRelativeTolerance = 1e-10;
final double parRelativeTolerance = 1e-10;
final double orthoTolerance = 1e-10;
final double threshold = Precision.SAFE_MIN;
// Extract the parameters to be fitted
final double[] initialSolution = getInitialSolution(a);
// TODO - Pass in more advanced stopping criteria.
// Create the target and weight arrays
final double[] yd = new double[n];
final double[] w = new double[n];
for (int i = 0; i < n; i++) {
yd[i] = y[i];
w[i] = 1;
}
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
//@formatter:off
LeastSquaresBuilder builder = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(getMaxEvaluations()).start(initialSolution).target(yd).weight(new DiagonalMatrix(w));
if (f instanceof ExtendedNonLinearFunction && ((ExtendedNonLinearFunction) f).canComputeValuesAndJacobian()) {
// Compute together, or each individually
builder.model(new ValueAndJacobianFunction() {
final ExtendedNonLinearFunction fun = (ExtendedNonLinearFunction) f;
public Pair<RealVector, RealMatrix> value(RealVector point) {
final double[] p = point.toArray();
final Pair<double[], double[][]> result = fun.computeValuesAndJacobian(p);
return new Pair<RealVector, RealMatrix>(new ArrayRealVector(result.getFirst(), false), new Array2DRowRealMatrix(result.getSecond(), false));
}
public RealVector computeValue(double[] params) {
return new ArrayRealVector(fun.computeValues(params), false);
}
public RealMatrix computeJacobian(double[] params) {
return new Array2DRowRealMatrix(fun.computeJacobian(params), false);
}
});
} else {
// Compute separately
builder.model(new MultivariateVectorFunctionWrapper((NonLinearFunction) f, a, n), new MultivariateMatrixFunctionWrapper((NonLinearFunction) f, a, n));
}
LeastSquaresProblem problem = builder.build();
Optimum optimum = optimizer.optimize(problem);
final double[] parameters = optimum.getPoint().toArray();
setSolution(a, parameters);
iterations = optimum.getIterations();
evaluations = optimum.getEvaluations();
if (a_dev != null) {
try {
double[][] covar = optimum.getCovariances(threshold).getData();
setDeviationsFromMatrix(a_dev, covar);
} catch (SingularMatrixException e) {
// Matrix inversion failed. In order to return a solution
// return the reciprocal of the diagonal of the Fisher information
// for a loose bound on the limit
final int[] gradientIndices = f.gradientIndices();
final int nparams = gradientIndices.length;
GradientCalculator calculator = GradientCalculatorFactory.newCalculator(nparams);
double[][] alpha = new double[nparams][nparams];
double[] beta = new double[nparams];
calculator.findLinearised(nparams, y, a, alpha, beta, (NonLinearFunction) f);
FisherInformationMatrix m = new FisherInformationMatrix(alpha);
setDeviations(a_dev, m.crlb(true));
}
}
// Compute function value
if (y_fit != null) {
Gaussian2DFunction f = (Gaussian2DFunction) this.f;
f.initialise0(a);
f.forEach(new ValueProcedure() {
int i = 0;
public void execute(double value) {
y_fit[i] = value;
}
});
}
// As this is unweighted then we can do this to get the sum of squared residuals
// This is the same as optimum.getCost() * optimum.getCost(); The getCost() function
// just computes the dot product anyway.
value = optimum.getResiduals().dotProduct(optimum.getResiduals());
} catch (TooManyEvaluationsException e) {
return FitStatus.TOO_MANY_EVALUATIONS;
} catch (TooManyIterationsException e) {
return FitStatus.TOO_MANY_ITERATIONS;
} catch (ConvergenceException e) {
// Occurs when QR decomposition fails - mark as a singular non-linear model (no solution)
return FitStatus.SINGULAR_NON_LINEAR_MODEL;
} catch (Exception e) {
// TODO - Find out the other exceptions from the fitter and add return values to match.
return FitStatus.UNKNOWN;
}
return FitStatus.OK;
}
use of org.apache.commons.math3.linear.SingularMatrixException in project knime-core by knime.
the class AbstractSGOptimizer method optimize.
public LogRegLearnerResult optimize(final int maxEpoch, final TrainingData<T> data, final Progress progress) throws CanceledExecutionException {
final int nRows = data.getRowCount();
final int nFets = data.getFeatureCount();
final int nCats = data.getTargetDimension();
final U updater = m_updaterFactory.create();
final WeightMatrix<T> beta = new SimpleWeightMatrix<>(nFets, nCats, true);
int epoch = 0;
for (; epoch < maxEpoch; epoch++) {
// notify learning rate strategy that a new epoch starts
m_lrStrategy.startNewEpoch(epoch);
progress.setProgress(((double) epoch) / maxEpoch, "Start epoch " + epoch + " of " + maxEpoch);
for (int k = 0; k < nRows; k++) {
progress.checkCanceled();
T x = data.getRandomRow();
prepareIteration(beta, x, updater, m_regUpdater, k);
double[] prediction = beta.predict(x);
double[] sig = m_loss.gradient(x, prediction);
double stepSize = m_lrStrategy.getCurrentLearningRate(x, prediction, sig);
// beta is updated in two steps
m_regUpdater.update(beta, stepSize, k);
performUpdate(x, updater, sig, beta, stepSize, k);
double scale = beta.getScale();
if (scale > 1e10 || scale < -1e10 || (scale > 0 && scale < 1e-10) || (scale < 0 && scale > -1e-10)) {
normalize(beta, updater, k);
beta.normalize();
}
}
postProcessEpoch(beta, updater, m_regUpdater);
if (m_stoppingCriterion.checkConvergence(beta)) {
break;
}
}
StringBuilder warnBuilder = new StringBuilder();
if (epoch >= maxEpoch) {
warnBuilder.append("The algorithm did not reach convergence after the specified number of epochs. " + "Setting the epoch limit higher might result in a better model.");
}
double lossSum = totalLoss(beta);
RealMatrix betaMat = MatrixUtils.createRealMatrix(beta.getWeightVector());
RealMatrix covMat = null;
if (m_calcCovMatrix) {
try {
covMat = calculateCovariateMatrix(beta);
} catch (SingularMatrixException e) {
if (warnBuilder.length() > 0) {
warnBuilder.append("\n");
}
warnBuilder.append("The covariance matrix could not be calculated because the" + " observed fisher information matrix was singular. Did you properly normalize the numerical features?");
covMat = null;
}
}
m_warning = warnBuilder.length() > 0 ? warnBuilder.toString() : null;
// in a maximum likelihood sense
return new LogRegLearnerResult(betaMat, covMat, epoch, -lossSum);
}
use of org.apache.commons.math3.linear.SingularMatrixException in project GDSC-SMLM by aherbert.
the class TrackPopulationAnalysis method createMixed.
/**
* Creates the multivariate gaussian mixture as the best of many repeats of the expectation
* maximisation algorithm.
*
* @param data the data
* @param dimensions the dimensions
* @param numComponents the number of components
* @return the multivariate gaussian mixture expectation maximization
*/
private MultivariateGaussianMixtureExpectationMaximization createMixed(final double[][] data, int dimensions, int numComponents) {
// Fit a mixed multivariate Gaussian with different repeats.
final UnitSphereSampler sampler = UnitSphereSampler.of(UniformRandomProviders.create(Mixers.stafford13(settings.seed++)), dimensions);
final LocalList<CompletableFuture<MultivariateGaussianMixtureExpectationMaximization>> results = new LocalList<>(settings.repeats);
final DoubleDoubleBiPredicate test = createConvergenceTest(settings.relativeError);
if (settings.debug) {
ImageJUtils.log(" Fitting %d components", numComponents);
}
final Ticker ticker = ImageJUtils.createTicker(settings.repeats, 2, "Fitting...");
final AtomicInteger failures = new AtomicInteger();
for (int i = 0; i < settings.repeats; i++) {
final double[] vector = sampler.sample();
results.add(CompletableFuture.supplyAsync(() -> {
final MultivariateGaussianMixtureExpectationMaximization fitter = new MultivariateGaussianMixtureExpectationMaximization(data);
try {
// This may also throw the same exceptions due to inversion of the covariance matrix
final MixtureMultivariateGaussianDistribution initialMixture = MultivariateGaussianMixtureExpectationMaximization.estimate(data, numComponents, point -> {
double dot = 0;
for (int j = 0; j < dimensions; j++) {
dot += vector[j] * point[j];
}
return dot;
});
final boolean result = fitter.fit(initialMixture, settings.maxIterations, test);
// Log the result. Note: The ImageJ log is synchronized.
if (settings.debug) {
ImageJUtils.log(" Fit: log-likelihood=%s, iter=%d, converged=%b", fitter.getLogLikelihood(), fitter.getIterations(), result);
}
return result ? fitter : null;
} catch (NonPositiveDefiniteMatrixException | SingularMatrixException ex) {
failures.getAndIncrement();
if (settings.debug) {
ImageJUtils.log(" Fit failed during iteration %d. No variance in a sub-population " + "component (check alpha is not always 1.0).", fitter.getIterations());
}
} finally {
ticker.tick();
}
return null;
}));
}
ImageJUtils.finished();
if (failures.get() != 0 && settings.debug) {
ImageJUtils.log(" %d component fit failed %d/%d", numComponents, failures.get(), settings.repeats);
}
// Collect results and return the best model.
return results.stream().map(f -> f.join()).filter(f -> f != null).sorted((f1, f2) -> Double.compare(f2.getLogLikelihood(), f1.getLogLikelihood())).findFirst().orElse(null);
}
use of org.apache.commons.math3.linear.SingularMatrixException in project gatk-protected by broadinstitute.
the class CoverageDropoutDetector method retrieveGaussianMixtureModelForFilteredTargets.
/** <p>Produces a Gaussian mixture model based on the difference between targets and segment means.</p>
* <p>Filters targets to populations where more than the minProportion lie in a single segment.</p>
* <p>Returns null if no pass filtering. Please note that in these cases,
* in the rest of this class, we use this to assume that a GMM is not a good model.</p>
*
* @param segments -- segments with segment mean in log2 copy ratio space
* @param targets -- targets with a log2 copy ratio estimate
* @param minProportion -- minimum proportion of all targets that a given segment must have in order to be used
* in the evaluation
* @param numComponents -- number of components to use in the GMM. Usually, this is 2.
* @return never {@code null}. Fitting result with indications whether it converged or was even attempted.
*/
private MixtureMultivariateNormalFitResult retrieveGaussianMixtureModelForFilteredTargets(final List<ModeledSegment> segments, final TargetCollection<ReadCountRecord.SingleSampleRecord> targets, double minProportion, int numComponents) {
// For each target in a segment that contains enough targets, normalize the difference against the segment mean
// and collapse the filtered targets into the copy ratio estimates.
final List<Double> filteredTargetsSegDiff = getNumProbeFilteredTargetList(segments, targets, minProportion);
if (filteredTargetsSegDiff.size() < numComponents) {
return new MixtureMultivariateNormalFitResult(null, false, false);
}
// Assume that Apache Commons wants data points in the first dimension.
// Note that second dimension of length 2 (instead of 1) is to wrok around funny Apache commons API.
final double[][] filteredTargetsSegDiff2d = new double[filteredTargetsSegDiff.size()][2];
// Convert the filtered targets into 2d array (even if second dimension is length 1). The second dimension is
// uncorrelated Gaussian. This is only to get around funny API in Apache Commons, which will throw an
// exception if the length of the second dimension is < 2
final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));
final NormalDistribution nd = new NormalDistribution(rng, 0, .1);
for (int i = 0; i < filteredTargetsSegDiff.size(); i++) {
filteredTargetsSegDiff2d[i][0] = filteredTargetsSegDiff.get(i);
filteredTargetsSegDiff2d[i][1] = nd.sample();
}
final MixtureMultivariateNormalDistribution estimateEM0 = MultivariateNormalMixtureExpectationMaximization.estimate(filteredTargetsSegDiff2d, numComponents);
final MultivariateNormalMixtureExpectationMaximization multivariateNormalMixtureExpectationMaximization = new MultivariateNormalMixtureExpectationMaximization(filteredTargetsSegDiff2d);
try {
multivariateNormalMixtureExpectationMaximization.fit(estimateEM0);
} catch (final MaxCountExceededException | ConvergenceException | SingularMatrixException e) {
// did not converge. Include the model as it was when the exception was thrown.
return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), false, true);
}
return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), true, true);
}
use of org.apache.commons.math3.linear.SingularMatrixException in project gatk by broadinstitute.
the class CoverageDropoutDetector method retrieveGaussianMixtureModelForFilteredTargets.
/** <p>Produces a Gaussian mixture model based on the difference between targets and segment means.</p>
* <p>Filters targets to populations where more than the minProportion lie in a single segment.</p>
* <p>Returns null if no pass filtering. Please note that in these cases,
* in the rest of this class, we use this to assume that a GMM is not a good model.</p>
*
* @param segments -- segments with segment mean in log2 copy ratio space
* @param targets -- targets with a log2 copy ratio estimate
* @param minProportion -- minimum proportion of all targets that a given segment must have in order to be used
* in the evaluation
* @param numComponents -- number of components to use in the GMM. Usually, this is 2.
* @return never {@code null}. Fitting result with indications whether it converged or was even attempted.
*/
private MixtureMultivariateNormalFitResult retrieveGaussianMixtureModelForFilteredTargets(final List<ModeledSegment> segments, final TargetCollection<ReadCountRecord.SingleSampleRecord> targets, double minProportion, int numComponents) {
// For each target in a segment that contains enough targets, normalize the difference against the segment mean
// and collapse the filtered targets into the copy ratio estimates.
final List<Double> filteredTargetsSegDiff = getNumProbeFilteredTargetList(segments, targets, minProportion);
if (filteredTargetsSegDiff.size() < numComponents) {
return new MixtureMultivariateNormalFitResult(null, false, false);
}
// Assume that Apache Commons wants data points in the first dimension.
// Note that second dimension of length 2 (instead of 1) is to wrok around funny Apache commons API.
final double[][] filteredTargetsSegDiff2d = new double[filteredTargetsSegDiff.size()][2];
// Convert the filtered targets into 2d array (even if second dimension is length 1). The second dimension is
// uncorrelated Gaussian. This is only to get around funny API in Apache Commons, which will throw an
// exception if the length of the second dimension is < 2
final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(RANDOM_SEED));
final NormalDistribution nd = new NormalDistribution(rng, 0, .1);
for (int i = 0; i < filteredTargetsSegDiff.size(); i++) {
filteredTargetsSegDiff2d[i][0] = filteredTargetsSegDiff.get(i);
filteredTargetsSegDiff2d[i][1] = nd.sample();
}
final MixtureMultivariateNormalDistribution estimateEM0 = MultivariateNormalMixtureExpectationMaximization.estimate(filteredTargetsSegDiff2d, numComponents);
final MultivariateNormalMixtureExpectationMaximization multivariateNormalMixtureExpectationMaximization = new MultivariateNormalMixtureExpectationMaximization(filteredTargetsSegDiff2d);
try {
multivariateNormalMixtureExpectationMaximization.fit(estimateEM0);
} catch (final MaxCountExceededException | ConvergenceException | SingularMatrixException e) {
// did not converge. Include the model as it was when the exception was thrown.
return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), false, true);
}
return new MixtureMultivariateNormalFitResult(multivariateNormalMixtureExpectationMaximization.getFittedModel(), true, true);
}
Aggregations