Implicit Weights
The weights are no longer implicit in LeastSquaresProblem.Evaluation. They are already included in the computed residuals and Jacobian. GN and LM multiplied the residuals by the weights immediately, so that was easy to remove. Created an AbstractEvaluation class which handles the derived quantitied (cost, rms, covariance,...) and two implementations. UnweightedEvaluation uses the straight forward formulas. DenseWeightedEvaluation delegates to an Evaluation and multiples the residuals and Jacobian by the square root of the weight matrix before returning them. Allowed me to remove the reference to the full weight matrix. git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1569344 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe
parent
567127c109
commit
06d490a4bd
|
|
@ -0,0 +1,70 @@
|
|||
package org.apache.commons.math3.fitting.leastsquares;
|
||||
|
||||
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation;
|
||||
import org.apache.commons.math3.linear.ArrayRealVector;
|
||||
import org.apache.commons.math3.linear.DecompositionSolver;
|
||||
import org.apache.commons.math3.linear.QRDecomposition;
|
||||
import org.apache.commons.math3.linear.RealMatrix;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
|
||||
/**
|
||||
* An implementation of {@link Evaluation} that is designed for extension. All of the
|
||||
* methods implemented here use the methods that are left unimplemented.
|
||||
* <p/>
|
||||
* TODO cache results?
|
||||
*
|
||||
* @version $Id$
|
||||
*/
|
||||
abstract class AbstractEvaluation implements Evaluation {
|
||||
|
||||
/** number of observations */
|
||||
private final int observationSize;
|
||||
|
||||
/**
|
||||
* Constructor.
|
||||
*
|
||||
* @param observationSize the number of observation. Needed for {@link
|
||||
* #computeRMS()}.
|
||||
*/
|
||||
AbstractEvaluation(final int observationSize) {
|
||||
this.observationSize = observationSize;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[][] computeCovariances(double threshold) {
|
||||
// Set up the Jacobian.
|
||||
final RealMatrix j = this.computeJacobian();
|
||||
|
||||
// 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();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] computeSigma(double covarianceSingularityThreshold) {
|
||||
final double[][] cov = this.computeCovariances(covarianceSingularityThreshold);
|
||||
final int nC = cov.length;
|
||||
final double[] sig = new double[nC];
|
||||
for (int i = 0; i < nC; ++i) {
|
||||
sig[i] = FastMath.sqrt(cov[i][i]);
|
||||
}
|
||||
return sig;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double computeRMS() {
|
||||
final double cost = this.computeCost();
|
||||
return FastMath.sqrt(cost * cost / this.observationSize);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double computeCost() {
|
||||
final ArrayRealVector r = new ArrayRealVector(this.computeResiduals());
|
||||
return FastMath.sqrt(r.dotProduct(r));
|
||||
}
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,55 @@
|
|||
package org.apache.commons.math3.fitting.leastsquares;
|
||||
|
||||
import org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation;
|
||||
import org.apache.commons.math3.linear.RealMatrix;
|
||||
|
||||
/**
|
||||
* Applies a dense weight matrix to an evaluation.
|
||||
*
|
||||
* @version $Id$
|
||||
*/
|
||||
class DenseWeightedEvaluation extends AbstractEvaluation {
|
||||
|
||||
/** the unweighted evaluation */
|
||||
private final Evaluation unweighted;
|
||||
/** reference to the weight square root matrix */
|
||||
private final RealMatrix weightSqrt;
|
||||
|
||||
/**
|
||||
* Create a weighted evaluation from an unweighted one.
|
||||
*
|
||||
* @param unweighted the evalutation before weights are applied
|
||||
* @param weightSqrt the matrix square root of the weight matrix
|
||||
*/
|
||||
DenseWeightedEvaluation(final Evaluation unweighted,
|
||||
final RealMatrix weightSqrt) {
|
||||
// weight square root is square, nR=nC=number of observations
|
||||
super(weightSqrt.getColumnDimension());
|
||||
this.unweighted = unweighted;
|
||||
this.weightSqrt = weightSqrt;
|
||||
}
|
||||
|
||||
/* apply weights */
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix computeJacobian() {
|
||||
return weightSqrt.multiply(this.unweighted.computeJacobian());
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] computeResiduals() {
|
||||
return this.weightSqrt.operate(this.unweighted.computeResiduals());
|
||||
}
|
||||
|
||||
/* delegate */
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] getPoint() {
|
||||
return unweighted.getPoint();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] computeValue() {
|
||||
return unweighted.computeValue();
|
||||
}
|
||||
}
|
||||
|
|
@ -104,17 +104,10 @@ public class GaussNewtonOptimizer implements LeastSquaresOptimizer {
|
|||
throw new NullArgumentException();
|
||||
}
|
||||
|
||||
final RealMatrix weightMatrix = lsp.getWeight();
|
||||
final int nR = weightMatrix.getRowDimension(); // Number of observed data.
|
||||
|
||||
// Diagonal of the weight matrix.
|
||||
final double[] residualsWeights = new double[nR];
|
||||
for (int i = 0; i < nR; i++) {
|
||||
residualsWeights[i] = weightMatrix.getEntry(i, i);
|
||||
}
|
||||
final int nR = lsp.getObservationSize(); // Number of observed data.
|
||||
final int nC = lsp.getParameterSize();
|
||||
|
||||
final double[] currentPoint = lsp.getStart();
|
||||
final int nC = currentPoint.length;
|
||||
|
||||
// iterate until convergence is reached
|
||||
PointVectorValuePair current = null;
|
||||
|
|
@ -128,7 +121,7 @@ public class GaussNewtonOptimizer implements LeastSquaresOptimizer {
|
|||
final Evaluation value = lsp.evaluate(currentPoint);
|
||||
final double[] currentObjective = value.computeValue();
|
||||
final double[] currentResiduals = value.computeResiduals();
|
||||
final RealMatrix weightedJacobian = value.computeWeightedJacobian();
|
||||
final RealMatrix weightedJacobian = value.computeJacobian();
|
||||
current = new PointVectorValuePair(currentPoint, currentObjective);
|
||||
|
||||
// build the linear problem
|
||||
|
|
@ -137,21 +130,20 @@ public class GaussNewtonOptimizer implements LeastSquaresOptimizer {
|
|||
for (int i = 0; i < nR; ++i) {
|
||||
|
||||
final double[] grad = weightedJacobian.getRow(i);
|
||||
final double weight = residualsWeights[i];
|
||||
final double residual = currentResiduals[i];
|
||||
|
||||
// compute the normal equation
|
||||
final double wr = weight * residual;
|
||||
//residual is already weighted
|
||||
for (int j = 0; j < nC; ++j) {
|
||||
b[j] += wr * grad[j];
|
||||
b[j] += residual * grad[j];
|
||||
}
|
||||
|
||||
// build the contribution matrix for measurement i
|
||||
for (int k = 0; k < nC; ++k) {
|
||||
double[] ak = a[k];
|
||||
double wgk = weight * grad[k];
|
||||
//Jacobian/gradient is already weighted
|
||||
for (int l = 0; l < nC; ++l) {
|
||||
ak[l] += wgk * grad[l];
|
||||
ak[l] += grad[k] * grad[l];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,60 @@
|
|||
package org.apache.commons.math3.fitting.leastsquares;
|
||||
|
||||
import org.apache.commons.math3.optim.ConvergenceChecker;
|
||||
import org.apache.commons.math3.optim.PointVectorValuePair;
|
||||
import org.apache.commons.math3.util.Incrementor;
|
||||
|
||||
/**
|
||||
* An adapter that delegates to another implementation of {@link LeastSquaresProblem}.
|
||||
*
|
||||
* @version $Id$
|
||||
*/
|
||||
public class LeastSquaresAdapter implements LeastSquaresProblem {
|
||||
|
||||
/** the delegate problem */
|
||||
private final LeastSquaresProblem problem;
|
||||
|
||||
/**
|
||||
* Delegate the {@link LeastSquaresProblem} interface to the given implementation.
|
||||
*
|
||||
* @param problem the delegate
|
||||
*/
|
||||
public LeastSquaresAdapter(final LeastSquaresProblem problem) {
|
||||
this.problem = problem;
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[] getStart() {
|
||||
return problem.getStart();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public int getObservationSize() {
|
||||
return problem.getObservationSize();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public int getParameterSize() {
|
||||
return problem.getParameterSize();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public Evaluation evaluate(final double[] point) {
|
||||
return problem.evaluate(point);
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public Incrementor getEvaluationCounter() {
|
||||
return problem.getEvaluationCounter();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public Incrementor getIterationCounter() {
|
||||
return problem.getIterationCounter();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public ConvergenceChecker<PointVectorValuePair> getConvergenceChecker() {
|
||||
return problem.getConvergenceChecker();
|
||||
}
|
||||
}
|
||||
|
|
@ -2,14 +2,57 @@ package org.apache.commons.math3.fitting.leastsquares;
|
|||
|
||||
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
|
||||
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
|
||||
import org.apache.commons.math3.linear.DiagonalMatrix;
|
||||
import org.apache.commons.math3.linear.EigenDecomposition;
|
||||
import org.apache.commons.math3.linear.RealMatrix;
|
||||
import org.apache.commons.math3.linear.RealVector;
|
||||
import org.apache.commons.math3.optim.ConvergenceChecker;
|
||||
import org.apache.commons.math3.optim.PointVectorValuePair;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
import org.apache.commons.math3.util.Incrementor;
|
||||
|
||||
/** @author Evan Ward */
|
||||
/**
|
||||
* A Factory for creating {@link LeastSquaresProblem}s.
|
||||
*
|
||||
* @version $Id$
|
||||
*/
|
||||
public class LeastSquaresFactory {
|
||||
|
||||
/** Prevent instantiation. */
|
||||
private LeastSquaresFactory() {
|
||||
}
|
||||
|
||||
/**
|
||||
* Create a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem}
|
||||
* from the given elements. There will be no weights applied (Identity weights).
|
||||
*
|
||||
* @param model the model function. Produces the computed values.
|
||||
* @param jacobian the jacobian of the model with respect to the parameters
|
||||
* @param observed the observed (target) values
|
||||
* @param start the initial guess.
|
||||
* @param checker convergence checker
|
||||
* @param maxEvaluations the maximum number of times to evaluate the model
|
||||
* @param maxIterations the maximum number to times to iterate in the algorithm
|
||||
* @return the specified General Least Squares problem.
|
||||
*/
|
||||
public static LeastSquaresProblem create(final MultivariateVectorFunction model,
|
||||
final MultivariateMatrixFunction jacobian,
|
||||
final double[] observed,
|
||||
final double[] start,
|
||||
final ConvergenceChecker<PointVectorValuePair> checker,
|
||||
final int maxEvaluations,
|
||||
final int maxIterations) {
|
||||
return new LeastSquaresProblemImpl(
|
||||
maxEvaluations,
|
||||
maxIterations,
|
||||
checker,
|
||||
observed,
|
||||
model,
|
||||
jacobian,
|
||||
start
|
||||
);
|
||||
}
|
||||
|
||||
/**
|
||||
* Create a {@link org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem}
|
||||
* from the given elements.
|
||||
|
|
@ -32,16 +75,50 @@ public class LeastSquaresFactory {
|
|||
final ConvergenceChecker<PointVectorValuePair> checker,
|
||||
final int maxEvaluations,
|
||||
final int maxIterations) {
|
||||
return new LeastSquaresProblemImpl(
|
||||
maxEvaluations,
|
||||
maxIterations,
|
||||
checker,
|
||||
observed,
|
||||
weight,
|
||||
model,
|
||||
jacobian,
|
||||
start
|
||||
);
|
||||
return weightMatrix(
|
||||
create(
|
||||
model,
|
||||
jacobian,
|
||||
observed,
|
||||
start,
|
||||
checker,
|
||||
maxEvaluations,
|
||||
maxIterations
|
||||
),
|
||||
weight);
|
||||
}
|
||||
|
||||
/**
|
||||
* Apply a dense weight matrix to the {@link LeastSquaresProblem}.
|
||||
*
|
||||
* @param problem the unweighted problem
|
||||
* @param weights the matrix of weights
|
||||
* @return a new {@link LeastSquaresProblem} with the weights applied. The original
|
||||
* {@code problem} is not modified.
|
||||
*/
|
||||
public static LeastSquaresProblem weightMatrix(final LeastSquaresProblem problem,
|
||||
final RealMatrix weights) {
|
||||
final RealMatrix weightSquareRoot = squareRoot(weights);
|
||||
return new LeastSquaresAdapter(problem) {
|
||||
@Override
|
||||
public Evaluation evaluate(final double[] point) {
|
||||
return new DenseWeightedEvaluation(super.evaluate(point), weightSquareRoot);
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
/**
|
||||
* Apply a diagon weight matrix to the {@link LeastSquaresProblem}.
|
||||
*
|
||||
* @param problem the unweighted problem
|
||||
* @param weights the diagonal of the weight matrix
|
||||
* @return a new {@link LeastSquaresProblem} with the weights applied. The original
|
||||
* {@code problem} is not modified.
|
||||
*/
|
||||
public static LeastSquaresProblem weightDiagonal(final LeastSquaresProblem problem,
|
||||
final RealVector weights) {
|
||||
//TODO more efficient implementation
|
||||
return weightMatrix(problem, new DiagonalMatrix(weights.toArray()));
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -55,48 +132,36 @@ public class LeastSquaresFactory {
|
|||
*/
|
||||
public static LeastSquaresProblem countEvaluations(final LeastSquaresProblem problem,
|
||||
final Incrementor counter) {
|
||||
//TODO adapter?
|
||||
return new LeastSquaresProblem() {
|
||||
return new LeastSquaresAdapter(problem) {
|
||||
|
||||
public Evaluation evaluate(double[] point) {
|
||||
public Evaluation evaluate(final double[] point) {
|
||||
counter.incrementCount();
|
||||
return problem.evaluate(point);
|
||||
return super.evaluate(point);
|
||||
}
|
||||
|
||||
/* delegate the rest */
|
||||
|
||||
public double[] getStart() {
|
||||
return problem.getStart();
|
||||
}
|
||||
|
||||
public int getObservationSize() {
|
||||
return problem.getObservationSize();
|
||||
}
|
||||
|
||||
public int getParameterSize() {
|
||||
return problem.getParameterSize();
|
||||
}
|
||||
|
||||
public RealMatrix getWeight() {
|
||||
return problem.getWeight();
|
||||
}
|
||||
|
||||
public RealMatrix getWeightSquareRoot() {
|
||||
return problem.getWeightSquareRoot();
|
||||
}
|
||||
|
||||
public Incrementor getEvaluationCounter() {
|
||||
return problem.getEvaluationCounter();
|
||||
}
|
||||
|
||||
public Incrementor getIterationCounter() {
|
||||
return problem.getIterationCounter();
|
||||
}
|
||||
|
||||
public ConvergenceChecker<PointVectorValuePair> getConvergenceChecker() {
|
||||
return problem.getConvergenceChecker();
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
/**
|
||||
* 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();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -39,22 +39,6 @@ public interface LeastSquaresProblem extends OptimizationProblem<PointVectorValu
|
|||
*/
|
||||
Evaluation evaluate(double[] point);
|
||||
|
||||
/**
|
||||
* Gets the weight matrix of the observations.
|
||||
* <p/>
|
||||
* TODO Is it possible to leave out this method and have the weights implicit in the
|
||||
* {@link Evaluation}?
|
||||
*
|
||||
* @return the weight matrix.
|
||||
*/
|
||||
RealMatrix getWeight();
|
||||
|
||||
/**Get the square root of the weight matrix.
|
||||
* TODO delete this method
|
||||
* @return the square root of the weight matrix
|
||||
*/
|
||||
RealMatrix getWeightSquareRoot();
|
||||
|
||||
public interface Evaluation {
|
||||
|
||||
/**
|
||||
|
|
@ -111,14 +95,7 @@ public interface LeastSquaresProblem extends OptimizationProblem<PointVectorValu
|
|||
* @throws DimensionMismatchException if the Jacobian dimension does not match
|
||||
* problem dimension.
|
||||
*/
|
||||
RealMatrix computeWeightedJacobian();
|
||||
|
||||
/**
|
||||
* Computes the Jacobian matrix.
|
||||
*
|
||||
* @return the Jacobian at the specified point.
|
||||
*/
|
||||
double[][] computeJacobian();
|
||||
RealMatrix computeJacobian();
|
||||
|
||||
/**
|
||||
* Computes the cost.
|
||||
|
|
@ -129,16 +106,16 @@ public interface LeastSquaresProblem extends OptimizationProblem<PointVectorValu
|
|||
double computeCost();
|
||||
|
||||
/**
|
||||
* Computes the residuals. The residual is the difference between the observed
|
||||
* (target) values and the model (objective function) value. There is one residual
|
||||
* for each element of the vector-valued function.
|
||||
* Computes the weighted residuals. The residual is the difference between the
|
||||
* observed (target) values and the model (objective function) value. There is one
|
||||
* residual for each element of the vector-valued function. The raw residuals are
|
||||
* then multiplied by the square root of the weight matrix.
|
||||
*
|
||||
* @return the residuals.
|
||||
* @return the weighted residuals: W<sup>1/2</sup> K.
|
||||
* @throws DimensionMismatchException if {@code params} has a wrong length.
|
||||
*/
|
||||
double[] computeResiduals();
|
||||
|
||||
|
||||
/**
|
||||
* Get the abscissa (independent variables) of this evaluation.
|
||||
*
|
||||
|
|
|
|||
|
|
@ -19,17 +19,11 @@ package org.apache.commons.math3.fitting.leastsquares;
|
|||
import org.apache.commons.math3.analysis.MultivariateMatrixFunction;
|
||||
import org.apache.commons.math3.analysis.MultivariateVectorFunction;
|
||||
import org.apache.commons.math3.exception.DimensionMismatchException;
|
||||
import org.apache.commons.math3.linear.ArrayRealVector;
|
||||
import org.apache.commons.math3.linear.DecompositionSolver;
|
||||
import org.apache.commons.math3.linear.DiagonalMatrix;
|
||||
import org.apache.commons.math3.linear.EigenDecomposition;
|
||||
import org.apache.commons.math3.linear.MatrixUtils;
|
||||
import org.apache.commons.math3.linear.QRDecomposition;
|
||||
import org.apache.commons.math3.linear.RealMatrix;
|
||||
import org.apache.commons.math3.optim.AbstractOptimizationProblem;
|
||||
import org.apache.commons.math3.optim.ConvergenceChecker;
|
||||
import org.apache.commons.math3.optim.PointVectorValuePair;
|
||||
import org.apache.commons.math3.util.FastMath;
|
||||
|
||||
/**
|
||||
* A private, "field" immutable (not "real" immutable) implementation of {@link
|
||||
|
|
@ -44,14 +38,10 @@ class LeastSquaresProblemImpl
|
|||
|
||||
/** Target values for the model function at optimum. */
|
||||
private double[] target;
|
||||
/** Weight matrix. */
|
||||
private RealMatrix weight;
|
||||
/** Model function. */
|
||||
private MultivariateVectorFunction model;
|
||||
/** Jacobian of the model function. */
|
||||
private MultivariateMatrixFunction jacobian;
|
||||
/** Square-root of the weight matrix. */
|
||||
private RealMatrix weightSqrt;
|
||||
/** Initial guess. */
|
||||
private double[] start;
|
||||
|
||||
|
|
@ -59,16 +49,13 @@ class LeastSquaresProblemImpl
|
|||
final int maxIterations,
|
||||
final ConvergenceChecker<PointVectorValuePair> checker,
|
||||
final double[] target,
|
||||
final RealMatrix weight,
|
||||
final MultivariateVectorFunction model,
|
||||
final MultivariateMatrixFunction jacobian,
|
||||
final double[] start) {
|
||||
super(maxEvaluations, maxIterations, checker);
|
||||
this.target = target;
|
||||
this.weight = weight;
|
||||
this.model = model;
|
||||
this.jacobian = jacobian;
|
||||
this.weightSqrt = squareRoot(weight);
|
||||
this.start = start;
|
||||
}
|
||||
|
||||
|
|
@ -93,44 +80,11 @@ class LeastSquaresProblemImpl
|
|||
return start == null ? null : start.clone();
|
||||
}
|
||||
|
||||
/**
|
||||
* Gets the square-root of the weight matrix.
|
||||
*
|
||||
* @return the square-root of the weight matrix.
|
||||
*/
|
||||
public RealMatrix getWeightSquareRoot() {
|
||||
return weightSqrt == null ? null : weightSqrt.copy();
|
||||
}
|
||||
|
||||
/**
|
||||
* Gets the model function.
|
||||
*
|
||||
* @return the model function.
|
||||
*/
|
||||
public MultivariateVectorFunction getModel() {
|
||||
return model;
|
||||
}
|
||||
|
||||
/**
|
||||
* Gets the model function's Jacobian.
|
||||
*
|
||||
* @return the Jacobian.
|
||||
*/
|
||||
public MultivariateMatrixFunction getJacobian() {
|
||||
return jacobian;
|
||||
}
|
||||
|
||||
public RealMatrix getWeight() {
|
||||
return weight.copy();
|
||||
}
|
||||
|
||||
public Evaluation evaluate(final double[] point) {
|
||||
//TODO evaluate value and jacobian in one function call
|
||||
return new EvaluationImpl(
|
||||
return new UnweightedEvaluation(
|
||||
this.model.value(point),
|
||||
this.jacobian.value(point),
|
||||
this.weight,
|
||||
this.weightSqrt,
|
||||
this.target,
|
||||
point);
|
||||
}
|
||||
|
|
@ -140,7 +94,7 @@ class LeastSquaresProblemImpl
|
|||
* <p/>
|
||||
* TODO revisit lazy evaluation
|
||||
*/
|
||||
private static class EvaluationImpl implements Evaluation {
|
||||
private static class UnweightedEvaluation extends AbstractEvaluation {
|
||||
|
||||
/** the point of evaluation */
|
||||
private final double[] point;
|
||||
|
|
@ -148,70 +102,32 @@ class LeastSquaresProblemImpl
|
|||
private final double[] values;
|
||||
/** deriviative at point */
|
||||
private final double[][] jacobian;
|
||||
/* references to data defined by the least squares problem. Not modified.
|
||||
* Could be a reference to the problem.
|
||||
*/
|
||||
private final RealMatrix weight;
|
||||
private final RealMatrix weightSqrt;
|
||||
/** reference to the observed values */
|
||||
private final double[] target;
|
||||
|
||||
private EvaluationImpl(final double[] values,
|
||||
final double[][] jacobian,
|
||||
final RealMatrix weight,
|
||||
final RealMatrix weightSqrt,
|
||||
final double[] target,
|
||||
final double[] point) {
|
||||
private UnweightedEvaluation(final double[] values,
|
||||
final double[][] jacobian,
|
||||
final double[] target,
|
||||
final double[] point) {
|
||||
super(target.length);
|
||||
this.values = values;
|
||||
this.jacobian = jacobian;
|
||||
this.weight = weight;
|
||||
this.weightSqrt = weightSqrt;
|
||||
this.target = target;
|
||||
this.point = point;
|
||||
}
|
||||
|
||||
public double[][] computeCovariances(double threshold) {
|
||||
// Set up the Jacobian.
|
||||
final RealMatrix j = computeWeightedJacobian();
|
||||
|
||||
// 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();
|
||||
}
|
||||
|
||||
public double[] computeSigma(double covarianceSingularityThreshold) {
|
||||
final double[][] cov = computeCovariances(covarianceSingularityThreshold);
|
||||
final int nC = cov.length;
|
||||
final double[] sig = new double[nC];
|
||||
for (int i = 0; i < nC; ++i) {
|
||||
sig[i] = FastMath.sqrt(cov[i][i]);
|
||||
}
|
||||
return sig;
|
||||
}
|
||||
|
||||
public double computeRMS() {
|
||||
final double cost = computeCost();
|
||||
return FastMath.sqrt(cost * cost / target.length);
|
||||
}
|
||||
|
||||
public double[] computeValue() {
|
||||
return this.values;
|
||||
}
|
||||
|
||||
public RealMatrix computeWeightedJacobian() {
|
||||
return weightSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian()));
|
||||
public RealMatrix computeJacobian() {
|
||||
return MatrixUtils.createRealMatrix(this.jacobian);
|
||||
}
|
||||
|
||||
public double[][] computeJacobian() {
|
||||
return this.jacobian;
|
||||
}
|
||||
|
||||
public double computeCost() {
|
||||
final ArrayRealVector r = new ArrayRealVector(computeResiduals());
|
||||
return FastMath.sqrt(r.dotProduct(weight.operate(r)));
|
||||
public double[] getPoint() {
|
||||
return this.point;
|
||||
}
|
||||
|
||||
public double[] computeResiduals() {
|
||||
|
|
@ -229,29 +145,6 @@ class LeastSquaresProblemImpl
|
|||
return residuals;
|
||||
}
|
||||
|
||||
public double[] getPoint() {
|
||||
//TODO copy?
|
||||
return this.point;
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Computes the square-root of the weight matrix.
|
||||
*
|
||||
* @param m Symmetric, positive-definite (weight) matrix.
|
||||
* @return the square-root of the weight matrix.
|
||||
*/
|
||||
private RealMatrix squareRoot(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();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -311,7 +311,6 @@ public class LevenbergMarquardtOptimizer implements LeastSquaresOptimizer {
|
|||
//convergence criterion
|
||||
final ConvergenceChecker<PointVectorValuePair> checker
|
||||
= problem.getConvergenceChecker();
|
||||
final RealMatrix weightMatrixSqrt = problem.getWeightSquareRoot();
|
||||
|
||||
// arrays shared with the other private methods
|
||||
final int solvedCols = FastMath.min(nR, nC);
|
||||
|
|
@ -350,13 +349,14 @@ public class LevenbergMarquardtOptimizer implements LeastSquaresOptimizer {
|
|||
|
||||
// QR decomposition of the jacobian matrix
|
||||
final InternalData internalData
|
||||
= qrDecomposition(value.computeWeightedJacobian(), solvedCols);
|
||||
= qrDecomposition(value.computeJacobian(), solvedCols);
|
||||
final double[][] weightedJacobian = internalData.weightedJacobian;
|
||||
final int[] permutation = internalData.permutation;
|
||||
final double[] diagR = internalData.diagR;
|
||||
final double[] jacNorm = internalData.jacNorm;
|
||||
|
||||
double[] weightedResidual = weightMatrixSqrt.operate(currentResiduals);
|
||||
//residuals already have weights applied
|
||||
double[] weightedResidual = currentResiduals;
|
||||
for (int i = 0; i < nR; i++) {
|
||||
qtf[i] = weightedResidual[i];
|
||||
}
|
||||
|
|
|
|||
|
|
@ -64,12 +64,7 @@ class OptimumImpl implements Optimum {
|
|||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public RealMatrix computeWeightedJacobian() {
|
||||
return value.computeWeightedJacobian();
|
||||
}
|
||||
|
||||
/** {@inheritDoc} */
|
||||
public double[][] computeJacobian() {
|
||||
public RealMatrix computeJacobian() {
|
||||
return value.computeJacobian();
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue