001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.fitting.leastsquares;
018
019import org.apache.commons.math4.exception.ConvergenceException;
020import org.apache.commons.math4.exception.NullArgumentException;
021import org.apache.commons.math4.exception.util.LocalizedFormats;
022import org.apache.commons.math4.fitting.leastsquares.LeastSquaresProblem.Evaluation;
023import org.apache.commons.math4.linear.ArrayRealVector;
024import org.apache.commons.math4.linear.CholeskyDecomposition;
025import org.apache.commons.math4.linear.LUDecomposition;
026import org.apache.commons.math4.linear.MatrixUtils;
027import org.apache.commons.math4.linear.NonPositiveDefiniteMatrixException;
028import org.apache.commons.math4.linear.QRDecomposition;
029import org.apache.commons.math4.linear.RealMatrix;
030import org.apache.commons.math4.linear.RealVector;
031import org.apache.commons.math4.linear.SingularMatrixException;
032import org.apache.commons.math4.linear.SingularValueDecomposition;
033import org.apache.commons.math4.optim.ConvergenceChecker;
034import org.apache.commons.math4.util.IntegerSequence;
035import org.apache.commons.math4.util.Pair;
036
037/**
038 * Gauss-Newton least-squares solver.
039 * <p> This class solve a least-square problem by
040 * solving the normal equations of the linearized problem at each iteration. Either LU
041 * decomposition or Cholesky decomposition can be used to solve the normal equations,
042 * or QR decomposition or SVD decomposition can be used to solve the linear system. LU
043 * decomposition is faster but QR decomposition is more robust for difficult problems,
044 * and SVD can compute a solution for rank-deficient problems.
045 * </p>
046 *
047 * @since 3.3
048 */
049public class GaussNewtonOptimizer implements LeastSquaresOptimizer {
050
051    /** The decomposition algorithm to use to solve the normal equations. */
052    //TODO move to linear package and expand options?
053    public enum Decomposition {
054        /**
055         * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and
056         * using the {@link LUDecomposition}.
057         *
058         * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the
059         * normal matrix and n<sup>3</sup>/3 operations (m &gt; n) to solve the system using
060         * the LU decomposition. </p>
061         */
062        LU {
063            @Override
064            protected RealVector solve(final RealMatrix jacobian,
065                                       final RealVector residuals) {
066                try {
067                    final Pair<RealMatrix, RealVector> normalEquation =
068                            computeNormalMatrix(jacobian, residuals);
069                    final RealMatrix normal = normalEquation.getFirst();
070                    final RealVector jTr = normalEquation.getSecond();
071                    return new LUDecomposition(normal, SINGULARITY_THRESHOLD)
072                            .getSolver()
073                            .solve(jTr);
074                } catch (SingularMatrixException e) {
075                    throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
076                }
077            }
078        },
079        /**
080         * Solve the linear least squares problem (Jx=r) using the {@link
081         * QRDecomposition}.
082         *
083         * <p> Theoretically this method takes mn<sup>2</sup> - n<sup>3</sup>/3 operations
084         * (m &gt; n) and has better numerical accuracy than any method that forms the normal
085         * equations. </p>
086         */
087        QR {
088            @Override
089            protected RealVector solve(final RealMatrix jacobian,
090                                       final RealVector residuals) {
091                try {
092                    return new QRDecomposition(jacobian, SINGULARITY_THRESHOLD)
093                            .getSolver()
094                            .solve(residuals);
095                } catch (SingularMatrixException e) {
096                    throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
097                }
098            }
099        },
100        /**
101         * Solve by forming the normal equations (J<sup>T</sup>Jx=J<sup>T</sup>r) and
102         * using the {@link CholeskyDecomposition}.
103         *
104         * <p> Theoretically this method takes mn<sup>2</sup>/2 operations to compute the
105         * normal matrix and n<sup>3</sup>/6 operations (m &gt; n) to solve the system using
106         * the Cholesky decomposition. </p>
107         */
108        CHOLESKY {
109            @Override
110            protected RealVector solve(final RealMatrix jacobian,
111                                       final RealVector residuals) {
112                try {
113                    final Pair<RealMatrix, RealVector> normalEquation =
114                            computeNormalMatrix(jacobian, residuals);
115                    final RealMatrix normal = normalEquation.getFirst();
116                    final RealVector jTr = normalEquation.getSecond();
117                    return new CholeskyDecomposition(
118                            normal, SINGULARITY_THRESHOLD, SINGULARITY_THRESHOLD)
119                            .getSolver()
120                            .solve(jTr);
121                } catch (NonPositiveDefiniteMatrixException e) {
122                    throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
123                }
124            }
125        },
126        /**
127         * Solve the linear least squares problem using the {@link
128         * SingularValueDecomposition}.
129         *
130         * <p> This method is slower, but can provide a solution for rank deficient and
131         * nearly singular systems.
132         */
133        SVD {
134            @Override
135            protected RealVector solve(final RealMatrix jacobian,
136                                       final RealVector residuals) {
137                return new SingularValueDecomposition(jacobian)
138                        .getSolver()
139                        .solve(residuals);
140            }
141        };
142
143        /**
144         * Solve the linear least squares problem Jx=r.
145         *
146         * @param jacobian  the Jacobian matrix, J. the number of rows &gt;= the number or
147         *                  columns.
148         * @param residuals the computed residuals, r.
149         * @return the solution x, to the linear least squares problem Jx=r.
150         * @throws ConvergenceException if the matrix properties (e.g. singular) do not
151         *                              permit a solution.
152         */
153        protected abstract RealVector solve(RealMatrix jacobian,
154                                            RealVector residuals);
155    }
156
157    /**
158     * The singularity threshold for matrix decompositions. Determines when a {@link
159     * ConvergenceException} is thrown. The current value was the default value for {@link
160     * LUDecomposition}.
161     */
162    private static final double SINGULARITY_THRESHOLD = 1e-11;
163
164    /** Indicator for using LU decomposition. */
165    private final Decomposition decomposition;
166
167    /**
168     * Creates a Gauss Newton optimizer.
169     * <p>
170     * The default for the algorithm is to solve the normal equations using QR
171     * decomposition.
172     */
173    public GaussNewtonOptimizer() {
174        this(Decomposition.QR);
175    }
176
177    /**
178     * Create a Gauss Newton optimizer that uses the given decomposition algorithm to
179     * solve the normal equations.
180     *
181     * @param decomposition the {@link Decomposition} algorithm.
182     */
183    public GaussNewtonOptimizer(final Decomposition decomposition) {
184        this.decomposition = decomposition;
185    }
186
187    /**
188     * Get the matrix decomposition algorithm used to solve the normal equations.
189     *
190     * @return the matrix {@link Decomposition} algoritm.
191     */
192    public Decomposition getDecomposition() {
193        return this.decomposition;
194    }
195
196    /**
197     * Configure the decomposition algorithm.
198     *
199     * @param newDecomposition the {@link Decomposition} algorithm to use.
200     * @return a new instance.
201     */
202    public GaussNewtonOptimizer withDecomposition(final Decomposition newDecomposition) {
203        return new GaussNewtonOptimizer(newDecomposition);
204    }
205
206    /** {@inheritDoc} */
207    @Override
208    public Optimum optimize(final LeastSquaresProblem lsp) {
209        //create local evaluation and iteration counts
210        final IntegerSequence.Incrementor evaluationCounter = lsp.getEvaluationCounter();
211        final IntegerSequence.Incrementor iterationCounter = lsp.getIterationCounter();
212        final ConvergenceChecker<Evaluation> checker
213                = lsp.getConvergenceChecker();
214
215        // Computation will be useless without a checker (see "for-loop").
216        if (checker == null) {
217            throw new NullArgumentException();
218        }
219
220        RealVector currentPoint = lsp.getStart();
221
222        // iterate until convergence is reached
223        Evaluation current = null;
224        while (true) {
225            iterationCounter.increment();
226
227            // evaluate the objective function and its jacobian
228            Evaluation previous = current;
229            // Value of the objective function at "currentPoint".
230            evaluationCounter.increment();
231            current = lsp.evaluate(currentPoint);
232            final RealVector currentResiduals = current.getResiduals();
233            final RealMatrix weightedJacobian = current.getJacobian();
234            currentPoint = current.getPoint();
235
236            // Check convergence.
237            if (previous != null &&
238                checker.converged(iterationCounter.getCount(), previous, current)) {
239                return new OptimumImpl(current,
240                                       evaluationCounter.getCount(),
241                                       iterationCounter.getCount());
242            }
243
244            // solve the linearized least squares problem
245            final RealVector dX = this.decomposition.solve(weightedJacobian, currentResiduals);
246            // update the estimated parameters
247            currentPoint = currentPoint.add(dX);
248        }
249    }
250
251    /** {@inheritDoc} */
252    @Override
253    public String toString() {
254        return "GaussNewtonOptimizer{" +
255                "decomposition=" + decomposition +
256                '}';
257    }
258
259    /**
260     * Compute the normal matrix, J<sup>T</sup>J.
261     *
262     * @param jacobian  the m by n jacobian matrix, J. Input.
263     * @param residuals the m by 1 residual vector, r. Input.
264     * @return  the n by n normal matrix and  the n by 1 J<sup>Tr vector.
265     */
266    private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian,
267                                                                    final RealVector residuals) {
268        //since the normal matrix is symmetric, we only need to compute half of it.
269        final int nR = jacobian.getRowDimension();
270        final int nC = jacobian.getColumnDimension();
271        //allocate space for return values
272        final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC);
273        final RealVector jTr = new ArrayRealVector(nC);
274        //for each measurement
275        for (int i = 0; i < nR; ++i) {
276            //compute JTr for measurement i
277            for (int j = 0; j < nC; j++) {
278                jTr.setEntry(j, jTr.getEntry(j) +
279                        residuals.getEntry(i) * jacobian.getEntry(i, j));
280            }
281
282            // add the contribution to the normal matrix for measurement i
283            for (int k = 0; k < nC; ++k) {
284                //only compute the upper triangular part
285                for (int l = k; l < nC; ++l) {
286                    normal.setEntry(k, l, normal.getEntry(k, l) +
287                            jacobian.getEntry(i, k) * jacobian.getEntry(i, l));
288                }
289            }
290        }
291        //copy the upper triangular part to the lower triangular part.
292        for (int i = 0; i < nC; i++) {
293            for (int j = 0; j < i; j++) {
294                normal.setEntry(i, j, normal.getEntry(j, i));
295            }
296        }
297        return new Pair<>(normal, jTr);
298    }
299
300}