http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java b/src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java deleted file mode 100644 index d64b442..0000000 --- a/src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java +++ /dev/null @@ -1,235 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.linear; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.MaxCountExceededException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.util.ExceptionContext; -import org.apache.commons.math3.util.IterationManager; - -/** - * <p> - * This is an implementation of the conjugate gradient method for - * {@link RealLinearOperator}. It follows closely the template by <a - * href="#BARR1994">Barrett et al. (1994)</a> (figure 2.5). The linear system at - * hand is A · x = b, and the residual is r = b - A · x. - * </p> - * <h3><a id="stopcrit">Default stopping criterion</a></h3> - * <p> - * A default stopping criterion is implemented. The iterations stop when || r || - * ≤ δ || b ||, where b is the right-hand side vector, r the current - * estimate of the residual, and δ a user-specified tolerance. It should - * be noted that r is the so-called <em>updated</em> residual, which might - * differ from the true residual due to rounding-off errors (see e.g. <a - * href="#STRA2002">Strakos and Tichy, 2002</a>). - * </p> - * <h3>Iteration count</h3> - * <p> - * In the present context, an iteration should be understood as one evaluation - * of the matrix-vector product A · x. The initialization phase therefore - * counts as one iteration. - * </p> - * <h3><a id="context">Exception context</a></h3> - * <p> - * Besides standard {@link DimensionMismatchException}, this class might throw - * {@link NonPositiveDefiniteOperatorException} if the linear operator or - * the preconditioner are not positive definite. In this case, the - * {@link ExceptionContext} provides some more information - * <ul> - * <li>key {@code "operator"} points to the offending linear operator, say L,</li> - * <li>key {@code "vector"} points to the offending vector, say x, such that - * x<sup>T</sup> · L · x < 0.</li> - * </ul> - * </p> - * <h3>References</h3> - * <dl> - * <dt><a id="BARR1994">Barret et al. (1994)</a></dt> - * <dd>R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra, - * V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst, - * <a href="http://www.netlib.org/linalg/html_templates/Templates.html";><em> - * Templates for the Solution of Linear Systems: Building Blocks for Iterative - * Methods</em></a>, SIAM</dd> - * <dt><a id="STRA2002">Strakos and Tichy (2002) - * <dt> - * <dd>Z. Strakos and P. Tichy, <a - * href="http://etna.mcs.kent.edu/vol.13.2002/pp56-80.dir/pp56-80.pdf";> - * <em>On error estimation in the conjugate gradient method and why it works - * in finite precision computations</em></a>, Electronic Transactions on - * Numerical Analysis 13: 56-80, 2002</dd> - * </dl> - * - * @since 3.0 - */ -public class ConjugateGradient - extends PreconditionedIterativeLinearSolver { - - /** Key for the <a href="#context">exception context</a>. */ - public static final String OPERATOR = "operator"; - - /** Key for the <a href="#context">exception context</a>. */ - public static final String VECTOR = "vector"; - - /** - * {@code true} if positive-definiteness of matrix and preconditioner should - * be checked. - */ - private boolean check; - - /** The value of δ, for the default stopping criterion. */ - private final double delta; - - /** - * Creates a new instance of this class, with <a href="#stopcrit">default - * stopping criterion</a>. - * - * @param maxIterations the maximum number of iterations - * @param delta the δ parameter for the default stopping criterion - * @param check {@code true} if positive definiteness of both matrix and - * preconditioner should be checked - */ - public ConjugateGradient(final int maxIterations, final double delta, - final boolean check) { - super(maxIterations); - this.delta = delta; - this.check = check; - } - - /** - * Creates a new instance of this class, with <a href="#stopcrit">default - * stopping criterion</a> and custom iteration manager. - * - * @param manager the custom iteration manager - * @param delta the δ parameter for the default stopping criterion - * @param check {@code true} if positive definiteness of both matrix and - * preconditioner should be checked - * @throws NullArgumentException if {@code manager} is {@code null} - */ - public ConjugateGradient(final IterationManager manager, - final double delta, final boolean check) - throws NullArgumentException { - super(manager); - this.delta = delta; - this.check = check; - } - - /** - * Returns {@code true} if positive-definiteness should be checked for both - * matrix and preconditioner. - * - * @return {@code true} if the tests are to be performed - */ - public final boolean getCheck() { - return check; - } - - /** - * {@inheritDoc} - * - * @throws NonPositiveDefiniteOperatorException if {@code a} or {@code m} is - * not positive definite - */ - @Override - public RealVector solveInPlace(final RealLinearOperator a, - final RealLinearOperator m, - final RealVector b, - final RealVector x0) - throws NullArgumentException, NonPositiveDefiniteOperatorException, - NonSquareOperatorException, DimensionMismatchException, - MaxCountExceededException { - checkParameters(a, m, b, x0); - final IterationManager manager = getIterationManager(); - // Initialization of default stopping criterion - manager.resetIterationCount(); - final double rmax = delta * b.getNorm(); - final RealVector bro = RealVector.unmodifiableRealVector(b); - - // Initialization phase counts as one iteration. - manager.incrementIterationCount(); - // p and x are constructed as copies of x0, since presumably, the type - // of x is optimized for the calculation of the matrix-vector product - // A.x. - final RealVector x = x0; - final RealVector xro = RealVector.unmodifiableRealVector(x); - final RealVector p = x.copy(); - RealVector q = a.operate(p); - - final RealVector r = b.combine(1, -1, q); - final RealVector rro = RealVector.unmodifiableRealVector(r); - double rnorm = r.getNorm(); - RealVector z; - if (m == null) { - z = r; - } else { - z = null; - } - IterativeLinearSolverEvent evt; - evt = new DefaultIterativeLinearSolverEvent(this, - manager.getIterations(), xro, bro, rro, rnorm); - manager.fireInitializationEvent(evt); - if (rnorm <= rmax) { - manager.fireTerminationEvent(evt); - return x; - } - double rhoPrev = 0.; - while (true) { - manager.incrementIterationCount(); - evt = new DefaultIterativeLinearSolverEvent(this, - manager.getIterations(), xro, bro, rro, rnorm); - manager.fireIterationStartedEvent(evt); - if (m != null) { - z = m.operate(r); - } - final double rhoNext = r.dotProduct(z); - if (check && (rhoNext <= 0.)) { - final NonPositiveDefiniteOperatorException e; - e = new NonPositiveDefiniteOperatorException(); - final ExceptionContext context = e.getContext(); - context.setValue(OPERATOR, m); - context.setValue(VECTOR, r); - throw e; - } - if (manager.getIterations() == 2) { - p.setSubVector(0, z); - } else { - p.combineToSelf(rhoNext / rhoPrev, 1., z); - } - q = a.operate(p); - final double pq = p.dotProduct(q); - if (check && (pq <= 0.)) { - final NonPositiveDefiniteOperatorException e; - e = new NonPositiveDefiniteOperatorException(); - final ExceptionContext context = e.getContext(); - context.setValue(OPERATOR, a); - context.setValue(VECTOR, p); - throw e; - } - final double alpha = rhoNext / pq; - x.combineToSelf(1., alpha, p); - r.combineToSelf(1., -alpha, q); - rhoPrev = rhoNext; - rnorm = r.getNorm(); - evt = new DefaultIterativeLinearSolverEvent(this, - manager.getIterations(), xro, bro, rro, rnorm); - manager.fireIterationPerformedEvent(evt); - if (rnorm <= rmax) { - manager.fireTerminationEvent(evt); - return x; - } - } - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/DecompositionSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/DecompositionSolver.java b/src/main/java/org/apache/commons/math3/linear/DecompositionSolver.java deleted file mode 100644 index 50090a9..0000000 --- a/src/main/java/org/apache/commons/math3/linear/DecompositionSolver.java +++ /dev/null @@ -1,97 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.linear; - -/** - * Interface handling decomposition algorithms that can solve A × X = B. - * <p> - * Decomposition algorithms decompose an A matrix has a product of several specific - * matrices from which they can solve A × X = B in least squares sense: they find X - * such that ||A × X - B|| is minimal. - * <p> - * Some solvers like {@link LUDecomposition} can only find the solution for - * square matrices and when the solution is an exact linear solution, i.e. when - * ||A × X - B|| is exactly 0. Other solvers can also find solutions - * with non-square matrix A and with non-null minimal norm. If an exact linear - * solution exists it is also the minimal norm solution. - * - * @since 2.0 - */ -public interface DecompositionSolver { - - /** - * Solve the linear equation A × X = B for matrices A. - * <p> - * The A matrix is implicit, it is provided by the underlying - * decomposition algorithm. - * - * @param b right-hand side of the equation A × X = B - * @return a vector X that minimizes the two norm of A × X - B - * @throws org.apache.commons.math3.exception.DimensionMismatchException - * if the matrices dimensions do not match. - * @throws SingularMatrixException if the decomposed matrix is singular. - */ - RealVector solve(final RealVector b) throws SingularMatrixException; - - /** - * Solve the linear equation A × X = B for matrices A. - * <p> - * The A matrix is implicit, it is provided by the underlying - * decomposition algorithm. - * - * @param b right-hand side of the equation A × X = B - * @return a matrix X that minimizes the two norm of A × X - B - * @throws org.apache.commons.math3.exception.DimensionMismatchException - * if the matrices dimensions do not match. - * @throws SingularMatrixException if the decomposed matrix is singular. - */ - RealMatrix solve(final RealMatrix b) throws SingularMatrixException; - - /** - * Check if the decomposed matrix is non-singular. - * @return true if the decomposed matrix is non-singular. - */ - boolean isNonSingular(); - - /** - * Get the <a href="http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse";>pseudo-inverse</a> - * of the decomposed matrix. - * <p> - * <em>This is equal to the inverse of the decomposed matrix, if such an inverse exists.</em> - * <p> - * If no such inverse exists, then the result has properties that resemble that of an inverse. - * <p> - * In particular, in this case, if the decomposed matrix is A, then the system of equations - * \( A x = b \) may have no solutions, or many. If it has no solutions, then the pseudo-inverse - * \( A^+ \) gives the "closest" solution \( z = A^+ b \), meaning \( \left \| A z - b \right \|_2 \) - * is minimized. If there are many solutions, then \( z = A^+ b \) is the smallest solution, - * meaning \( \left \| z \right \|_2 \) is minimized. - * <p> - * Note however that some decompositions cannot compute a pseudo-inverse for all matrices. - * For example, the {@link LUDecomposition} is not defined for non-square matrices to begin - * with. The {@link QRDecomposition} can operate on non-square matrices, but will throw - * {@link SingularMatrixException} if the decomposed matrix is singular. Refer to the javadoc - * of specific decomposition implementations for more details. - * - * @return pseudo-inverse matrix (which is the inverse, if it exists), - * if the decomposition can pseudo-invert the decomposed matrix - * @throws SingularMatrixException if the decomposed matrix is singular and the decomposition - * can not compute a pseudo-inverse - */ - RealMatrix getInverse() throws SingularMatrixException; -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/DefaultFieldMatrixChangingVisitor.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/DefaultFieldMatrixChangingVisitor.java b/src/main/java/org/apache/commons/math3/linear/DefaultFieldMatrixChangingVisitor.java deleted file mode 100644 index 11a61fa..0000000 --- a/src/main/java/org/apache/commons/math3/linear/DefaultFieldMatrixChangingVisitor.java +++ /dev/null @@ -1,58 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.linear; - -import org.apache.commons.math3.FieldElement; - -/** - * Default implementation of the {@link FieldMatrixChangingVisitor} interface. - * <p> - * This class is a convenience to create custom visitors without defining all - * methods. This class provides default implementations that do nothing. - * </p> - * - * @param <T> the type of the field elements - * @since 2.0 - */ -public class DefaultFieldMatrixChangingVisitor<T extends FieldElement<T>> - implements FieldMatrixChangingVisitor<T> { - /** Zero element of the field. */ - private final T zero; - - /** Build a new instance. - * @param zero additive identity of the field - */ - public DefaultFieldMatrixChangingVisitor(final T zero) { - this.zero = zero; - } - - /** {@inheritDoc} */ - public void start(int rows, int columns, - int startRow, int endRow, int startColumn, int endColumn) { - } - - /** {@inheritDoc} */ - public T visit(int row, int column, T value) { - return value; - } - - /** {@inheritDoc} */ - public T end() { - return zero; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/DefaultFieldMatrixPreservingVisitor.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/DefaultFieldMatrixPreservingVisitor.java b/src/main/java/org/apache/commons/math3/linear/DefaultFieldMatrixPreservingVisitor.java deleted file mode 100644 index 97fb59b..0000000 --- a/src/main/java/org/apache/commons/math3/linear/DefaultFieldMatrixPreservingVisitor.java +++ /dev/null @@ -1,56 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.linear; - -import org.apache.commons.math3.FieldElement; - -/** - * Default implementation of the {@link FieldMatrixPreservingVisitor} interface. - * <p> - * This class is a convenience to create custom visitors without defining all - * methods. This class provides default implementations that do nothing. - * </p> - * - * @param <T> the type of the field elements - * @since 2.0 - */ -public class DefaultFieldMatrixPreservingVisitor<T extends FieldElement<T>> - implements FieldMatrixPreservingVisitor<T> { - /** Zero element of the field. */ - private final T zero; - - /** Build a new instance. - * @param zero additive identity of the field - */ - public DefaultFieldMatrixPreservingVisitor(final T zero) { - this.zero = zero; - } - - /** {@inheritDoc} */ - public void start(int rows, int columns, - int startRow, int endRow, int startColumn, int endColumn) { - } - - /** {@inheritDoc} */ - public void visit(int row, int column, T value) {} - - /** {@inheritDoc} */ - public T end() { - return zero; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/DefaultIterativeLinearSolverEvent.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/DefaultIterativeLinearSolverEvent.java b/src/main/java/org/apache/commons/math3/linear/DefaultIterativeLinearSolverEvent.java deleted file mode 100644 index dbced15..0000000 --- a/src/main/java/org/apache/commons/math3/linear/DefaultIterativeLinearSolverEvent.java +++ /dev/null @@ -1,143 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.linear; - -import org.apache.commons.math3.exception.MathUnsupportedOperationException; - -/** - * A default concrete implementation of the abstract class - * {@link IterativeLinearSolverEvent}. - * - */ -public class DefaultIterativeLinearSolverEvent extends IterativeLinearSolverEvent { - - /** */ - private static final long serialVersionUID = 20120129L; - - /** The right-hand side vector. */ - private final RealVector b; - - /** The current estimate of the residual. */ - private final RealVector r; - - /** The current estimate of the norm of the residual. */ - private final double rnorm; - - /** The current estimate of the solution. */ - private final RealVector x; - - /** - * Creates a new instance of this class. This implementation does - * <em>not</em> deep copy the specified vectors {@code x}, {@code b}, - * {@code r}. Therefore the user must make sure that these vectors are - * either unmodifiable views or deep copies of the same vectors actually - * used by the {@code source}. Failure to do so may compromise subsequent - * iterations of the {@code source}. If the residual vector {@code r} is - * {@code null}, then {@link #getResidual()} throws a - * {@link MathUnsupportedOperationException}, and - * {@link #providesResidual()} returns {@code false}. - * - * @param source the iterative solver which fired this event - * @param iterations the number of iterations performed at the time - * {@code this} event is created - * @param x the current estimate of the solution - * @param b the right-hand side vector - * @param r the current estimate of the residual (can be {@code null}) - * @param rnorm the norm of the current estimate of the residual - */ - public DefaultIterativeLinearSolverEvent(final Object source, final int iterations, - final RealVector x, final RealVector b, final RealVector r, - final double rnorm) { - super(source, iterations); - this.x = x; - this.b = b; - this.r = r; - this.rnorm = rnorm; - } - - /** - * Creates a new instance of this class. This implementation does - * <em>not</em> deep copy the specified vectors {@code x}, {@code b}. - * Therefore the user must make sure that these vectors are either - * unmodifiable views or deep copies of the same vectors actually used by - * the {@code source}. Failure to do so may compromise subsequent iterations - * of the {@code source}. Callling {@link #getResidual()} on instances - * returned by this constructor throws a - * {@link MathUnsupportedOperationException}, while - * {@link #providesResidual()} returns {@code false}. - * - * @param source the iterative solver which fired this event - * @param iterations the number of iterations performed at the time - * {@code this} event is created - * @param x the current estimate of the solution - * @param b the right-hand side vector - * @param rnorm the norm of the current estimate of the residual - */ - public DefaultIterativeLinearSolverEvent(final Object source, final int iterations, - final RealVector x, final RealVector b, final double rnorm) { - super(source, iterations); - this.x = x; - this.b = b; - this.r = null; - this.rnorm = rnorm; - } - - /** {@inheritDoc} */ - @Override - public double getNormOfResidual() { - return rnorm; - } - - /** - * {@inheritDoc} - * - * This implementation throws an {@link MathUnsupportedOperationException} - * if no residual vector {@code r} was provided at construction time. - */ - @Override - public RealVector getResidual() { - if (r != null) { - return r; - } - throw new MathUnsupportedOperationException(); - } - - /** {@inheritDoc} */ - @Override - public RealVector getRightHandSideVector() { - return b; - } - - /** {@inheritDoc} */ - @Override - public RealVector getSolution() { - return x; - } - - /** - * {@inheritDoc} - * - * This implementation returns {@code true} if a non-{@code null} value was - * specified for the residual vector {@code r} at construction time. - * - * @return {@code true} if {@code r != null} - */ - @Override - public boolean providesResidual() { - return r != null; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/DefaultRealMatrixChangingVisitor.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/DefaultRealMatrixChangingVisitor.java b/src/main/java/org/apache/commons/math3/linear/DefaultRealMatrixChangingVisitor.java deleted file mode 100644 index 64949f7..0000000 --- a/src/main/java/org/apache/commons/math3/linear/DefaultRealMatrixChangingVisitor.java +++ /dev/null @@ -1,44 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.linear; - -/** - * Default implementation of the {@link RealMatrixChangingVisitor} interface. - * <p> - * This class is a convenience to create custom visitors without defining all - * methods. This class provides default implementations that do nothing. - * </p> - * - * @since 2.0 - */ -public class DefaultRealMatrixChangingVisitor implements RealMatrixChangingVisitor { - /** {@inheritDoc} */ - public void start(int rows, int columns, - int startRow, int endRow, int startColumn, int endColumn) { - } - - /** {@inheritDoc} */ - public double visit(int row, int column, double value) { - return value; - } - - /** {@inheritDoc} */ - public double end() { - return 0; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/DefaultRealMatrixPreservingVisitor.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/DefaultRealMatrixPreservingVisitor.java b/src/main/java/org/apache/commons/math3/linear/DefaultRealMatrixPreservingVisitor.java deleted file mode 100644 index 103e85b..0000000 --- a/src/main/java/org/apache/commons/math3/linear/DefaultRealMatrixPreservingVisitor.java +++ /dev/null @@ -1,42 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.linear; - -/** - * Default implementation of the {@link RealMatrixPreservingVisitor} interface. - * <p> - * This class is a convenience to create custom visitors without defining all - * methods. This class provides default implementations that do nothing. - * </p> - * - * @since 2.0 - */ -public class DefaultRealMatrixPreservingVisitor implements RealMatrixPreservingVisitor { - /** {@inheritDoc} */ - public void start(int rows, int columns, - int startRow, int endRow, int startColumn, int endColumn) { - } - - /** {@inheritDoc} */ - public void visit(int row, int column, double value) {} - - /** {@inheritDoc} */ - public double end() { - return 0; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/DiagonalMatrix.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/DiagonalMatrix.java b/src/main/java/org/apache/commons/math3/linear/DiagonalMatrix.java deleted file mode 100644 index 22ab9f1..0000000 --- a/src/main/java/org/apache/commons/math3/linear/DiagonalMatrix.java +++ /dev/null @@ -1,370 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.linear; - -import java.io.Serializable; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.NumberIsTooLargeException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.MathUtils; -import org.apache.commons.math3.util.Precision; - -/** - * Implementation of a diagonal matrix. - * - * @since 3.1.1 - */ -public class DiagonalMatrix extends AbstractRealMatrix - implements Serializable { - /** Serializable version identifier. */ - private static final long serialVersionUID = 20121229L; - /** Entries of the diagonal. */ - private final double[] data; - - /** - * Creates a matrix with the supplied dimension. - * - * @param dimension Number of rows and columns in the new matrix. - * @throws NotStrictlyPositiveException if the dimension is - * not positive. - */ - public DiagonalMatrix(final int dimension) - throws NotStrictlyPositiveException { - super(dimension, dimension); - data = new double[dimension]; - } - - /** - * Creates a matrix using the input array as the underlying data. - * <br/> - * The input array is copied, not referenced. - * - * @param d Data for the new matrix. - */ - public DiagonalMatrix(final double[] d) { - this(d, true); - } - - /** - * Creates a matrix using the input array as the underlying data. - * <br/> - * If an array is created specially in order to be embedded in a - * this instance and not used directly, the {@code copyArray} may be - * set to {@code false}. - * This will prevent the copying and improve performance as no new - * array will be built and no data will be copied. - * - * @param d Data for new matrix. - * @param copyArray if {@code true}, the input array will be copied, - * otherwise it will be referenced. - * @exception NullArgumentException if d is null - */ - public DiagonalMatrix(final double[] d, final boolean copyArray) - throws NullArgumentException { - MathUtils.checkNotNull(d); - data = copyArray ? d.clone() : d; - } - - /** - * {@inheritDoc} - * - * @throws DimensionMismatchException if the requested dimensions are not equal. - */ - @Override - public RealMatrix createMatrix(final int rowDimension, - final int columnDimension) - throws NotStrictlyPositiveException, - DimensionMismatchException { - if (rowDimension != columnDimension) { - throw new DimensionMismatchException(rowDimension, columnDimension); - } - - return new DiagonalMatrix(rowDimension); - } - - /** {@inheritDoc} */ - @Override - public RealMatrix copy() { - return new DiagonalMatrix(data); - } - - /** - * Compute the sum of {@code this} and {@code m}. - * - * @param m Matrix to be added. - * @return {@code this + m}. - * @throws MatrixDimensionMismatchException if {@code m} is not the same - * size as {@code this}. - */ - public DiagonalMatrix add(final DiagonalMatrix m) - throws MatrixDimensionMismatchException { - // Safety check. - MatrixUtils.checkAdditionCompatible(this, m); - - final int dim = getRowDimension(); - final double[] outData = new double[dim]; - for (int i = 0; i < dim; i++) { - outData[i] = data[i] + m.data[i]; - } - - return new DiagonalMatrix(outData, false); - } - - /** - * Returns {@code this} minus {@code m}. - * - * @param m Matrix to be subtracted. - * @return {@code this - m} - * @throws MatrixDimensionMismatchException if {@code m} is not the same - * size as {@code this}. - */ - public DiagonalMatrix subtract(final DiagonalMatrix m) - throws MatrixDimensionMismatchException { - MatrixUtils.checkSubtractionCompatible(this, m); - - final int dim = getRowDimension(); - final double[] outData = new double[dim]; - for (int i = 0; i < dim; i++) { - outData[i] = data[i] - m.data[i]; - } - - return new DiagonalMatrix(outData, false); - } - - /** - * Returns the result of postmultiplying {@code this} by {@code m}. - * - * @param m matrix to postmultiply by - * @return {@code this * m} - * @throws DimensionMismatchException if - * {@code columnDimension(this) != rowDimension(m)} - */ - public DiagonalMatrix multiply(final DiagonalMatrix m) - throws DimensionMismatchException { - MatrixUtils.checkMultiplicationCompatible(this, m); - - final int dim = getRowDimension(); - final double[] outData = new double[dim]; - for (int i = 0; i < dim; i++) { - outData[i] = data[i] * m.data[i]; - } - - return new DiagonalMatrix(outData, false); - } - - /** - * Returns the result of postmultiplying {@code this} by {@code m}. - * - * @param m matrix to postmultiply by - * @return {@code this * m} - * @throws DimensionMismatchException if - * {@code columnDimension(this) != rowDimension(m)} - */ - @Override - public RealMatrix multiply(final RealMatrix m) - throws DimensionMismatchException { - if (m instanceof DiagonalMatrix) { - return multiply((DiagonalMatrix) m); - } else { - MatrixUtils.checkMultiplicationCompatible(this, m); - final int nRows = m.getRowDimension(); - final int nCols = m.getColumnDimension(); - final double[][] product = new double[nRows][nCols]; - for (int r = 0; r < nRows; r++) { - for (int c = 0; c < nCols; c++) { - product[r][c] = data[r] * m.getEntry(r, c); - } - } - return new Array2DRowRealMatrix(product, false); - } - } - - /** {@inheritDoc} */ - @Override - public double[][] getData() { - final int dim = getRowDimension(); - final double[][] out = new double[dim][dim]; - - for (int i = 0; i < dim; i++) { - out[i][i] = data[i]; - } - - return out; - } - - /** - * Gets a reference to the underlying data array. - * - * @return 1-dimensional array of entries. - */ - public double[] getDataRef() { - return data; - } - - /** {@inheritDoc} */ - @Override - public double getEntry(final int row, final int column) - throws OutOfRangeException { - MatrixUtils.checkMatrixIndex(this, row, column); - return row == column ? data[row] : 0; - } - - /** {@inheritDoc} - * @throws NumberIsTooLargeException if {@code row != column} and value is non-zero. - */ - @Override - public void setEntry(final int row, final int column, final double value) - throws OutOfRangeException, NumberIsTooLargeException { - if (row == column) { - MatrixUtils.checkRowIndex(this, row); - data[row] = value; - } else { - ensureZero(value); - } - } - - /** {@inheritDoc} - * @throws NumberIsTooLargeException if {@code row != column} and increment is non-zero. - */ - @Override - public void addToEntry(final int row, - final int column, - final double increment) - throws OutOfRangeException, NumberIsTooLargeException { - if (row == column) { - MatrixUtils.checkRowIndex(this, row); - data[row] += increment; - } else { - ensureZero(increment); - } - } - - /** {@inheritDoc} */ - @Override - public void multiplyEntry(final int row, - final int column, - final double factor) - throws OutOfRangeException { - // we don't care about non-diagonal elements for multiplication - if (row == column) { - MatrixUtils.checkRowIndex(this, row); - data[row] *= factor; - } - } - - /** {@inheritDoc} */ - @Override - public int getRowDimension() { - return data.length; - } - - /** {@inheritDoc} */ - @Override - public int getColumnDimension() { - return data.length; - } - - /** {@inheritDoc} */ - @Override - public double[] operate(final double[] v) - throws DimensionMismatchException { - return multiply(new DiagonalMatrix(v, false)).getDataRef(); - } - - /** {@inheritDoc} */ - @Override - public double[] preMultiply(final double[] v) - throws DimensionMismatchException { - return operate(v); - } - - /** {@inheritDoc} */ - @Override - public RealVector preMultiply(final RealVector v) throws DimensionMismatchException { - final double[] vectorData; - if (v instanceof ArrayRealVector) { - vectorData = ((ArrayRealVector) v).getDataRef(); - } else { - vectorData = v.toArray(); - } - return MatrixUtils.createRealVector(preMultiply(vectorData)); - } - - /** Ensure a value is zero. - * @param value value to check - * @exception NumberIsTooLargeException if value is not zero - */ - private void ensureZero(final double value) throws NumberIsTooLargeException { - if (!Precision.equals(0.0, value, 1)) { - throw new NumberIsTooLargeException(FastMath.abs(value), 0, true); - } - } - - /** - * Computes the inverse of this diagonal matrix. - * <p> - * Note: this method will use a singularity threshold of 0, - * use {@link #inverse(double)} if a different threshold is needed. - * - * @return the inverse of {@code m} - * @throws SingularMatrixException if the matrix is singular - * @since 3.3 - */ - public DiagonalMatrix inverse() throws SingularMatrixException { - return inverse(0); - } - - /** - * Computes the inverse of this diagonal matrix. - * - * @param threshold Singularity threshold. - * @return the inverse of {@code m} - * @throws SingularMatrixException if the matrix is singular - * @since 3.3 - */ - public DiagonalMatrix inverse(double threshold) throws SingularMatrixException { - if (isSingular(threshold)) { - throw new SingularMatrixException(); - } - - final double[] result = new double[data.length]; - for (int i = 0; i < data.length; i++) { - result[i] = 1.0 / data[i]; - } - return new DiagonalMatrix(result, false); - } - - /** Returns whether this diagonal matrix is singular, i.e. any diagonal entry - * is equal to {@code 0} within the given threshold. - * - * @param threshold Singularity threshold. - * @return {@code true} if the matrix is singular, {@code false} otherwise - * @since 3.3 - */ - public boolean isSingular(double threshold) { - for (int i = 0; i < data.length; i++) { - if (Precision.equals(data[i], 0.0, threshold)) { - return true; - } - } - return false; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java b/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java deleted file mode 100644 index bd3819e..0000000 --- a/src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java +++ /dev/null @@ -1,970 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.linear; - -import org.apache.commons.math3.complex.Complex; -import org.apache.commons.math3.exception.MathArithmeticException; -import org.apache.commons.math3.exception.MathUnsupportedOperationException; -import org.apache.commons.math3.exception.MaxCountExceededException; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.util.Precision; -import org.apache.commons.math3.util.FastMath; - -/** - * Calculates the eigen decomposition of a real matrix. - * <p>The eigen decomposition of matrix A is a set of two matrices: - * V and D such that A = V × D × V<sup>T</sup>. - * A, V and D are all m × m matrices.</p> - * <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code> - * class from the <a href="http://math.nist.gov/javanumerics/jama/";>JAMA</a> - * library, with the following changes:</p> - * <ul> - * <li>a {@link #getVT() getVt} method has been added,</li> - * <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int) - * getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li> - * <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single - * eigenvector has been added,</li> - * <li>a {@link #getDeterminant() getDeterminant} method has been added.</li> - * <li>a {@link #getSolver() getSolver} method has been added.</li> - * </ul> - * <p> - * As of 3.1, this class supports general real matrices (both symmetric and non-symmetric): - * </p> - * <p> - * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal and the eigenvector - * matrix V is orthogonal, i.e. A = V.multiply(D.multiply(V.transpose())) and - * V.multiply(V.transpose()) equals the identity matrix. - * </p> - * <p> - * If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real eigenvalues - * in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2 blocks: - * <pre> - * [lambda, mu ] - * [ -mu, lambda] - * </pre> - * The columns of V represent the eigenvectors in the sense that A*V = V*D, - * i.e. A.multiply(V) equals V.multiply(D). - * The matrix V may be badly conditioned, or even singular, so the validity of the equation - * A = V*D*inverse(V) depends upon the condition of V. - * </p> - * <p> - * This implementation is based on the paper by A. Drubrulle, R.S. Martin and - * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971) - * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag, - * New-York - * </p> - * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html";>MathWorld</a> - * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix";>Wikipedia</a> - * @since 2.0 (changed to concrete class in 3.0) - */ -public class EigenDecomposition { - /** Internally used epsilon criteria. */ - private static final double EPSILON = 1e-12; - /** Maximum number of iterations accepted in the implicit QL transformation */ - private byte maxIter = 30; - /** Main diagonal of the tridiagonal matrix. */ - private double[] main; - /** Secondary diagonal of the tridiagonal matrix. */ - private double[] secondary; - /** - * Transformer to tridiagonal (may be null if matrix is already - * tridiagonal). - */ - private TriDiagonalTransformer transformer; - /** Real part of the realEigenvalues. */ - private double[] realEigenvalues; - /** Imaginary part of the realEigenvalues. */ - private double[] imagEigenvalues; - /** Eigenvectors. */ - private ArrayRealVector[] eigenvectors; - /** Cached value of V. */ - private RealMatrix cachedV; - /** Cached value of D. */ - private RealMatrix cachedD; - /** Cached value of Vt. */ - private RealMatrix cachedVt; - /** Whether the matrix is symmetric. */ - private final boolean isSymmetric; - - /** - * Calculates the eigen decomposition of the given real matrix. - * <p> - * Supports decomposition of a general matrix since 3.1. - * - * @param matrix Matrix to decompose. - * @throws MaxCountExceededException if the algorithm fails to converge. - * @throws MathArithmeticException if the decomposition of a general matrix - * results in a matrix with zero norm - * @since 3.1 - */ - public EigenDecomposition(final RealMatrix matrix) - throws MathArithmeticException { - final double symTol = 10 * matrix.getRowDimension() * matrix.getColumnDimension() * Precision.EPSILON; - isSymmetric = MatrixUtils.isSymmetric(matrix, symTol); - if (isSymmetric) { - transformToTridiagonal(matrix); - findEigenVectors(transformer.getQ().getData()); - } else { - final SchurTransformer t = transformToSchur(matrix); - findEigenVectorsFromSchur(t); - } - } - - /** - * Calculates the eigen decomposition of the given real matrix. - * - * @param matrix Matrix to decompose. - * @param splitTolerance Dummy parameter (present for backward - * compatibility only). - * @throws MathArithmeticException if the decomposition of a general matrix - * results in a matrix with zero norm - * @throws MaxCountExceededException if the algorithm fails to converge. - * @deprecated in 3.1 (to be removed in 4.0) due to unused parameter - */ - @Deprecated - public EigenDecomposition(final RealMatrix matrix, - final double splitTolerance) - throws MathArithmeticException { - this(matrix); - } - - /** - * Calculates the eigen decomposition of the symmetric tridiagonal - * matrix. The Householder matrix is assumed to be the identity matrix. - * - * @param main Main diagonal of the symmetric tridiagonal form. - * @param secondary Secondary of the tridiagonal form. - * @throws MaxCountExceededException if the algorithm fails to converge. - * @since 3.1 - */ - public EigenDecomposition(final double[] main, final double[] secondary) { - isSymmetric = true; - this.main = main.clone(); - this.secondary = secondary.clone(); - transformer = null; - final int size = main.length; - final double[][] z = new double[size][size]; - for (int i = 0; i < size; i++) { - z[i][i] = 1.0; - } - findEigenVectors(z); - } - - /** - * Calculates the eigen decomposition of the symmetric tridiagonal - * matrix. The Householder matrix is assumed to be the identity matrix. - * - * @param main Main diagonal of the symmetric tridiagonal form. - * @param secondary Secondary of the tridiagonal form. - * @param splitTolerance Dummy parameter (present for backward - * compatibility only). - * @throws MaxCountExceededException if the algorithm fails to converge. - * @deprecated in 3.1 (to be removed in 4.0) due to unused parameter - */ - @Deprecated - public EigenDecomposition(final double[] main, final double[] secondary, - final double splitTolerance) { - this(main, secondary); - } - - /** - * Gets the matrix V of the decomposition. - * V is an orthogonal matrix, i.e. its transpose is also its inverse. - * The columns of V are the eigenvectors of the original matrix. - * No assumption is made about the orientation of the system axes formed - * by the columns of V (e.g. in a 3-dimension space, V can form a left- - * or right-handed system). - * - * @return the V matrix. - */ - public RealMatrix getV() { - - if (cachedV == null) { - final int m = eigenvectors.length; - cachedV = MatrixUtils.createRealMatrix(m, m); - for (int k = 0; k < m; ++k) { - cachedV.setColumnVector(k, eigenvectors[k]); - } - } - // return the cached matrix - return cachedV; - } - - /** - * Gets the block diagonal matrix D of the decomposition. - * D is a block diagonal matrix. - * Real eigenvalues are on the diagonal while complex values are on - * 2x2 blocks { {real +imaginary}, {-imaginary, real} }. - * - * @return the D matrix. - * - * @see #getRealEigenvalues() - * @see #getImagEigenvalues() - */ - public RealMatrix getD() { - - if (cachedD == null) { - // cache the matrix for subsequent calls - cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues); - - for (int i = 0; i < imagEigenvalues.length; i++) { - if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) > 0) { - cachedD.setEntry(i, i+1, imagEigenvalues[i]); - } else if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) { - cachedD.setEntry(i, i-1, imagEigenvalues[i]); - } - } - } - return cachedD; - } - - /** - * Gets the transpose of the matrix V of the decomposition. - * V is an orthogonal matrix, i.e. its transpose is also its inverse. - * The columns of V are the eigenvectors of the original matrix. - * No assumption is made about the orientation of the system axes formed - * by the columns of V (e.g. in a 3-dimension space, V can form a left- - * or right-handed system). - * - * @return the transpose of the V matrix. - */ - public RealMatrix getVT() { - - if (cachedVt == null) { - final int m = eigenvectors.length; - cachedVt = MatrixUtils.createRealMatrix(m, m); - for (int k = 0; k < m; ++k) { - cachedVt.setRowVector(k, eigenvectors[k]); - } - } - - // return the cached matrix - return cachedVt; - } - - /** - * Returns whether the calculated eigen values are complex or real. - * <p>The method performs a zero check for each element of the - * {@link #getImagEigenvalues()} array and returns {@code true} if any - * element is not equal to zero. - * - * @return {@code true} if the eigen values are complex, {@code false} otherwise - * @since 3.1 - */ - public boolean hasComplexEigenvalues() { - for (int i = 0; i < imagEigenvalues.length; i++) { - if (!Precision.equals(imagEigenvalues[i], 0.0, EPSILON)) { - return true; - } - } - return false; - } - - /** - * Gets a copy of the real parts of the eigenvalues of the original matrix. - * - * @return a copy of the real parts of the eigenvalues of the original matrix. - * - * @see #getD() - * @see #getRealEigenvalue(int) - * @see #getImagEigenvalues() - */ - public double[] getRealEigenvalues() { - return realEigenvalues.clone(); - } - - /** - * Returns the real part of the i<sup>th</sup> eigenvalue of the original - * matrix. - * - * @param i index of the eigenvalue (counting from 0) - * @return real part of the i<sup>th</sup> eigenvalue of the original - * matrix. - * - * @see #getD() - * @see #getRealEigenvalues() - * @see #getImagEigenvalue(int) - */ - public double getRealEigenvalue(final int i) { - return realEigenvalues[i]; - } - - /** - * Gets a copy of the imaginary parts of the eigenvalues of the original - * matrix. - * - * @return a copy of the imaginary parts of the eigenvalues of the original - * matrix. - * - * @see #getD() - * @see #getImagEigenvalue(int) - * @see #getRealEigenvalues() - */ - public double[] getImagEigenvalues() { - return imagEigenvalues.clone(); - } - - /** - * Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original - * matrix. - * - * @param i Index of the eigenvalue (counting from 0). - * @return the imaginary part of the i<sup>th</sup> eigenvalue of the original - * matrix. - * - * @see #getD() - * @see #getImagEigenvalues() - * @see #getRealEigenvalue(int) - */ - public double getImagEigenvalue(final int i) { - return imagEigenvalues[i]; - } - - /** - * Gets a copy of the i<sup>th</sup> eigenvector of the original matrix. - * - * @param i Index of the eigenvector (counting from 0). - * @return a copy of the i<sup>th</sup> eigenvector of the original matrix. - * @see #getD() - */ - public RealVector getEigenvector(final int i) { - return eigenvectors[i].copy(); - } - - /** - * Computes the determinant of the matrix. - * - * @return the determinant of the matrix. - */ - public double getDeterminant() { - double determinant = 1; - for (double lambda : realEigenvalues) { - determinant *= lambda; - } - return determinant; - } - - /** - * Computes the square-root of the matrix. - * This implementation assumes that the matrix is symmetric and positive - * definite. - * - * @return the square-root of the matrix. - * @throws MathUnsupportedOperationException if the matrix is not - * symmetric or not positive definite. - * @since 3.1 - */ - public RealMatrix getSquareRoot() { - if (!isSymmetric) { - throw new MathUnsupportedOperationException(); - } - - final double[] sqrtEigenValues = new double[realEigenvalues.length]; - for (int i = 0; i < realEigenvalues.length; i++) { - final double eigen = realEigenvalues[i]; - if (eigen <= 0) { - throw new MathUnsupportedOperationException(); - } - sqrtEigenValues[i] = FastMath.sqrt(eigen); - } - final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues); - final RealMatrix v = getV(); - final RealMatrix vT = getVT(); - - return v.multiply(sqrtEigen).multiply(vT); - } - - /** - * Gets a solver for finding the A × X = B solution in exact - * linear sense. - * <p> - * Since 3.1, eigen decomposition of a general matrix is supported, - * but the {@link DecompositionSolver} only supports real eigenvalues. - * - * @return a solver - * @throws MathUnsupportedOperationException if the decomposition resulted in - * complex eigenvalues - */ - public DecompositionSolver getSolver() { - if (hasComplexEigenvalues()) { - throw new MathUnsupportedOperationException(); - } - return new Solver(realEigenvalues, imagEigenvalues, eigenvectors); - } - - /** Specialized solver. */ - private static class Solver implements DecompositionSolver { - /** Real part of the realEigenvalues. */ - private double[] realEigenvalues; - /** Imaginary part of the realEigenvalues. */ - private double[] imagEigenvalues; - /** Eigenvectors. */ - private final ArrayRealVector[] eigenvectors; - - /** - * Builds a solver from decomposed matrix. - * - * @param realEigenvalues Real parts of the eigenvalues. - * @param imagEigenvalues Imaginary parts of the eigenvalues. - * @param eigenvectors Eigenvectors. - */ - private Solver(final double[] realEigenvalues, - final double[] imagEigenvalues, - final ArrayRealVector[] eigenvectors) { - this.realEigenvalues = realEigenvalues; - this.imagEigenvalues = imagEigenvalues; - this.eigenvectors = eigenvectors; - } - - /** - * Solves the linear equation A × X = B for symmetric matrices A. - * <p> - * This method only finds exact linear solutions, i.e. solutions for - * which ||A × X - B|| is exactly 0. - * </p> - * - * @param b Right-hand side of the equation A × X = B. - * @return a Vector X that minimizes the two norm of A × X - B. - * - * @throws DimensionMismatchException if the matrices dimensions do not match. - * @throws SingularMatrixException if the decomposed matrix is singular. - */ - public RealVector solve(final RealVector b) { - if (!isNonSingular()) { - throw new SingularMatrixException(); - } - - final int m = realEigenvalues.length; - if (b.getDimension() != m) { - throw new DimensionMismatchException(b.getDimension(), m); - } - - final double[] bp = new double[m]; - for (int i = 0; i < m; ++i) { - final ArrayRealVector v = eigenvectors[i]; - final double[] vData = v.getDataRef(); - final double s = v.dotProduct(b) / realEigenvalues[i]; - for (int j = 0; j < m; ++j) { - bp[j] += s * vData[j]; - } - } - - return new ArrayRealVector(bp, false); - } - - /** {@inheritDoc} */ - public RealMatrix solve(RealMatrix b) { - - if (!isNonSingular()) { - throw new SingularMatrixException(); - } - - final int m = realEigenvalues.length; - 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; - } - for (int i = 0; i < m; ++i) { - final ArrayRealVector v = eigenvectors[i]; - final double[] vData = v.getDataRef(); - double s = 0; - for (int j = 0; j < m; ++j) { - s += v.getEntry(j) * tmpCol[j]; - } - s /= realEigenvalues[i]; - for (int j = 0; j < m; ++j) { - bp[j][k] += s * vData[j]; - } - } - } - - return new Array2DRowRealMatrix(bp, false); - - } - - /** - * Checks whether the decomposed matrix is non-singular. - * - * @return true if the decomposed matrix is non-singular. - */ - public boolean isNonSingular() { - double largestEigenvalueNorm = 0.0; - // Looping over all values (in case they are not sorted in decreasing - // order of their norm). - for (int i = 0; i < realEigenvalues.length; ++i) { - largestEigenvalueNorm = FastMath.max(largestEigenvalueNorm, eigenvalueNorm(i)); - } - // Corner case: zero matrix, all exactly 0 eigenvalues - if (largestEigenvalueNorm == 0.0) { - return false; - } - for (int i = 0; i < realEigenvalues.length; ++i) { - // Looking for eigenvalues that are 0, where we consider anything much much smaller - // than the largest eigenvalue to be effectively 0. - if (Precision.equals(eigenvalueNorm(i) / largestEigenvalueNorm, 0, EPSILON)) { - return false; - } - } - return true; - } - - /** - * @param i which eigenvalue to find the norm of - * @return the norm of ith (complex) eigenvalue. - */ - private double eigenvalueNorm(int i) { - final double re = realEigenvalues[i]; - final double im = imagEigenvalues[i]; - return FastMath.sqrt(re * re + im * im); - } - - /** - * Get the inverse of the decomposed matrix. - * - * @return the inverse matrix. - * @throws SingularMatrixException if the decomposed matrix is singular. - */ - public RealMatrix getInverse() { - if (!isNonSingular()) { - throw new SingularMatrixException(); - } - - final int m = realEigenvalues.length; - final double[][] invData = new double[m][m]; - - for (int i = 0; i < m; ++i) { - final double[] invI = invData[i]; - for (int j = 0; j < m; ++j) { - double invIJ = 0; - for (int k = 0; k < m; ++k) { - final double[] vK = eigenvectors[k].getDataRef(); - invIJ += vK[i] * vK[j] / realEigenvalues[k]; - } - invI[j] = invIJ; - } - } - return MatrixUtils.createRealMatrix(invData); - } - } - - /** - * Transforms the matrix to tridiagonal form. - * - * @param matrix Matrix to transform. - */ - private void transformToTridiagonal(final RealMatrix matrix) { - // transform the matrix to tridiagonal - transformer = new TriDiagonalTransformer(matrix); - main = transformer.getMainDiagonalRef(); - secondary = transformer.getSecondaryDiagonalRef(); - } - - /** - * Find eigenvalues and eigenvectors (Dubrulle et al., 1971) - * - * @param householderMatrix Householder matrix of the transformation - * to tridiagonal form. - */ - private void findEigenVectors(final double[][] householderMatrix) { - final double[][]z = householderMatrix.clone(); - final int n = main.length; - realEigenvalues = new double[n]; - imagEigenvalues = new double[n]; - final double[] e = new double[n]; - for (int i = 0; i < n - 1; i++) { - realEigenvalues[i] = main[i]; - e[i] = secondary[i]; - } - realEigenvalues[n - 1] = main[n - 1]; - e[n - 1] = 0; - - // Determine the largest main and secondary value in absolute term. - double maxAbsoluteValue = 0; - for (int i = 0; i < n; i++) { - if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) { - maxAbsoluteValue = FastMath.abs(realEigenvalues[i]); - } - if (FastMath.abs(e[i]) > maxAbsoluteValue) { - maxAbsoluteValue = FastMath.abs(e[i]); - } - } - // Make null any main and secondary value too small to be significant - if (maxAbsoluteValue != 0) { - for (int i=0; i < n; i++) { - if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) { - realEigenvalues[i] = 0; - } - if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) { - e[i]=0; - } - } - } - - for (int j = 0; j < n; j++) { - int its = 0; - int m; - do { - for (m = j; m < n - 1; m++) { - double delta = FastMath.abs(realEigenvalues[m]) + - FastMath.abs(realEigenvalues[m + 1]); - if (FastMath.abs(e[m]) + delta == delta) { - break; - } - } - if (m != j) { - if (its == maxIter) { - throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED, - maxIter); - } - its++; - double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]); - double t = FastMath.sqrt(1 + q * q); - if (q < 0.0) { - q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t); - } else { - q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t); - } - double u = 0.0; - double s = 1.0; - double c = 1.0; - int i; - for (i = m - 1; i >= j; i--) { - double p = s * e[i]; - double h = c * e[i]; - if (FastMath.abs(p) >= FastMath.abs(q)) { - c = q / p; - t = FastMath.sqrt(c * c + 1.0); - e[i + 1] = p * t; - s = 1.0 / t; - c *= s; - } else { - s = p / q; - t = FastMath.sqrt(s * s + 1.0); - e[i + 1] = q * t; - c = 1.0 / t; - s *= c; - } - if (e[i + 1] == 0.0) { - realEigenvalues[i + 1] -= u; - e[m] = 0.0; - break; - } - q = realEigenvalues[i + 1] - u; - t = (realEigenvalues[i] - q) * s + 2.0 * c * h; - u = s * t; - realEigenvalues[i + 1] = q + u; - q = c * t - h; - for (int ia = 0; ia < n; ia++) { - p = z[ia][i + 1]; - z[ia][i + 1] = s * z[ia][i] + c * p; - z[ia][i] = c * z[ia][i] - s * p; - } - } - if (t == 0.0 && i >= j) { - continue; - } - realEigenvalues[j] -= u; - e[j] = q; - e[m] = 0.0; - } - } while (m != j); - } - - //Sort the eigen values (and vectors) in increase order - for (int i = 0; i < n; i++) { - int k = i; - double p = realEigenvalues[i]; - for (int j = i + 1; j < n; j++) { - if (realEigenvalues[j] > p) { - k = j; - p = realEigenvalues[j]; - } - } - if (k != i) { - realEigenvalues[k] = realEigenvalues[i]; - realEigenvalues[i] = p; - for (int j = 0; j < n; j++) { - p = z[j][i]; - z[j][i] = z[j][k]; - z[j][k] = p; - } - } - } - - // Determine the largest eigen value in absolute term. - maxAbsoluteValue = 0; - for (int i = 0; i < n; i++) { - if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) { - maxAbsoluteValue=FastMath.abs(realEigenvalues[i]); - } - } - // Make null any eigen value too small to be significant - if (maxAbsoluteValue != 0.0) { - for (int i=0; i < n; i++) { - if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) { - realEigenvalues[i] = 0; - } - } - } - eigenvectors = new ArrayRealVector[n]; - final double[] tmp = new double[n]; - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - tmp[j] = z[j][i]; - } - eigenvectors[i] = new ArrayRealVector(tmp); - } - } - - /** - * Transforms the matrix to Schur form and calculates the eigenvalues. - * - * @param matrix Matrix to transform. - * @return the {@link SchurTransformer Shur transform} for this matrix - */ - private SchurTransformer transformToSchur(final RealMatrix matrix) { - final SchurTransformer schurTransform = new SchurTransformer(matrix); - final double[][] matT = schurTransform.getT().getData(); - - realEigenvalues = new double[matT.length]; - imagEigenvalues = new double[matT.length]; - - for (int i = 0; i < realEigenvalues.length; i++) { - if (i == (realEigenvalues.length - 1) || - Precision.equals(matT[i + 1][i], 0.0, EPSILON)) { - realEigenvalues[i] = matT[i][i]; - } else { - final double x = matT[i + 1][i + 1]; - final double p = 0.5 * (matT[i][i] - x); - final double z = FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1])); - realEigenvalues[i] = x + p; - imagEigenvalues[i] = z; - realEigenvalues[i + 1] = x + p; - imagEigenvalues[i + 1] = -z; - i++; - } - } - return schurTransform; - } - - /** - * Performs a division of two complex numbers. - * - * @param xr real part of the first number - * @param xi imaginary part of the first number - * @param yr real part of the second number - * @param yi imaginary part of the second number - * @return result of the complex division - */ - private Complex cdiv(final double xr, final double xi, - final double yr, final double yi) { - return new Complex(xr, xi).divide(new Complex(yr, yi)); - } - - /** - * Find eigenvectors from a matrix transformed to Schur form. - * - * @param schur the schur transformation of the matrix - * @throws MathArithmeticException if the Schur form has a norm of zero - */ - private void findEigenVectorsFromSchur(final SchurTransformer schur) - throws MathArithmeticException { - final double[][] matrixT = schur.getT().getData(); - final double[][] matrixP = schur.getP().getData(); - - final int n = matrixT.length; - - // compute matrix norm - double norm = 0.0; - for (int i = 0; i < n; i++) { - for (int j = FastMath.max(i - 1, 0); j < n; j++) { - norm += FastMath.abs(matrixT[i][j]); - } - } - - // we can not handle a matrix with zero norm - if (Precision.equals(norm, 0.0, EPSILON)) { - throw new MathArithmeticException(LocalizedFormats.ZERO_NORM); - } - - // Backsubstitute to find vectors of upper triangular form - - double r = 0.0; - double s = 0.0; - double z = 0.0; - - for (int idx = n - 1; idx >= 0; idx--) { - double p = realEigenvalues[idx]; - double q = imagEigenvalues[idx]; - - if (Precision.equals(q, 0.0)) { - // Real vector - int l = idx; - matrixT[idx][idx] = 1.0; - for (int i = idx - 1; i >= 0; i--) { - double w = matrixT[i][i] - p; - r = 0.0; - for (int j = l; j <= idx; j++) { - r += matrixT[i][j] * matrixT[j][idx]; - } - if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) { - z = w; - s = r; - } else { - l = i; - if (Precision.equals(imagEigenvalues[i], 0.0)) { - if (w != 0.0) { - matrixT[i][idx] = -r / w; - } else { - matrixT[i][idx] = -r / (Precision.EPSILON * norm); - } - } else { - // Solve real equations - double x = matrixT[i][i + 1]; - double y = matrixT[i + 1][i]; - q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) + - imagEigenvalues[i] * imagEigenvalues[i]; - double t = (x * s - z * r) / q; - matrixT[i][idx] = t; - if (FastMath.abs(x) > FastMath.abs(z)) { - matrixT[i + 1][idx] = (-r - w * t) / x; - } else { - matrixT[i + 1][idx] = (-s - y * t) / z; - } - } - - // Overflow control - double t = FastMath.abs(matrixT[i][idx]); - if ((Precision.EPSILON * t) * t > 1) { - for (int j = i; j <= idx; j++) { - matrixT[j][idx] /= t; - } - } - } - } - } else if (q < 0.0) { - // Complex vector - int l = idx - 1; - - // Last vector component imaginary so matrix is triangular - if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) { - matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1]; - matrixT[idx - 1][idx] = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1]; - } else { - final Complex result = cdiv(0.0, -matrixT[idx - 1][idx], - matrixT[idx - 1][idx - 1] - p, q); - matrixT[idx - 1][idx - 1] = result.getReal(); - matrixT[idx - 1][idx] = result.getImaginary(); - } - - matrixT[idx][idx - 1] = 0.0; - matrixT[idx][idx] = 1.0; - - for (int i = idx - 2; i >= 0; i--) { - double ra = 0.0; - double sa = 0.0; - for (int j = l; j <= idx; j++) { - ra += matrixT[i][j] * matrixT[j][idx - 1]; - sa += matrixT[i][j] * matrixT[j][idx]; - } - double w = matrixT[i][i] - p; - - if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) { - z = w; - r = ra; - s = sa; - } else { - l = i; - if (Precision.equals(imagEigenvalues[i], 0.0)) { - final Complex c = cdiv(-ra, -sa, w, q); - matrixT[i][idx - 1] = c.getReal(); - matrixT[i][idx] = c.getImaginary(); - } else { - // Solve complex equations - double x = matrixT[i][i + 1]; - double y = matrixT[i + 1][i]; - double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) + - imagEigenvalues[i] * imagEigenvalues[i] - q * q; - final double vi = (realEigenvalues[i] - p) * 2.0 * q; - if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) { - vr = Precision.EPSILON * norm * - (FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) + - FastMath.abs(y) + FastMath.abs(z)); - } - final Complex c = cdiv(x * r - z * ra + q * sa, - x * s - z * sa - q * ra, vr, vi); - matrixT[i][idx - 1] = c.getReal(); - matrixT[i][idx] = c.getImaginary(); - - if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) { - matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] + - q * matrixT[i][idx]) / x; - matrixT[i + 1][idx] = (-sa - w * matrixT[i][idx] - - q * matrixT[i][idx - 1]) / x; - } else { - final Complex c2 = cdiv(-r - y * matrixT[i][idx - 1], - -s - y * matrixT[i][idx], z, q); - matrixT[i + 1][idx - 1] = c2.getReal(); - matrixT[i + 1][idx] = c2.getImaginary(); - } - } - - // Overflow control - double t = FastMath.max(FastMath.abs(matrixT[i][idx - 1]), - FastMath.abs(matrixT[i][idx])); - if ((Precision.EPSILON * t) * t > 1) { - for (int j = i; j <= idx; j++) { - matrixT[j][idx - 1] /= t; - matrixT[j][idx] /= t; - } - } - } - } - } - } - - // Back transformation to get eigenvectors of original matrix - for (int j = n - 1; j >= 0; j--) { - for (int i = 0; i <= n - 1; i++) { - z = 0.0; - for (int k = 0; k <= FastMath.min(j, n - 1); k++) { - z += matrixP[i][k] * matrixT[k][j]; - } - matrixP[i][j] = z; - } - } - - eigenvectors = new ArrayRealVector[n]; - final double[] tmp = new double[n]; - for (int i = 0; i < n; i++) { - for (int j = 0; j < n; j++) { - tmp[j] = matrixP[j][i]; - } - eigenvectors[i] = new ArrayRealVector(tmp); - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/linear/FieldDecompositionSolver.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/linear/FieldDecompositionSolver.java b/src/main/java/org/apache/commons/math3/linear/FieldDecompositionSolver.java deleted file mode 100644 index 322eb3b..0000000 --- a/src/main/java/org/apache/commons/math3/linear/FieldDecompositionSolver.java +++ /dev/null @@ -1,75 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -package org.apache.commons.math3.linear; - -import org.apache.commons.math3.FieldElement; - - -/** - * Interface handling decomposition algorithms that can solve A × X = B. - * <p>Decomposition algorithms decompose an A matrix has a product of several specific - * matrices from which they can solve A × X = B in least squares sense: they find X - * such that ||A × X - B|| is minimal.</p> - * <p>Some solvers like {@link FieldLUDecomposition} can only find the solution for - * square matrices and when the solution is an exact linear solution, i.e. when - * ||A × X - B|| is exactly 0. Other solvers can also find solutions - * with non-square matrix A and with non-null minimal norm. If an exact linear - * solution exists it is also the minimal norm solution.</p> - * - * @param <T> the type of the field elements - * @since 2.0 - */ -public interface FieldDecompositionSolver<T extends FieldElement<T>> { - - /** Solve the linear equation A × X = B for matrices A. - * <p>The A matrix is implicit, it is provided by the underlying - * decomposition algorithm.</p> - * @param b right-hand side of the equation A × X = B - * @return a vector X that minimizes the two norm of A × X - B - * @throws org.apache.commons.math3.exception.DimensionMismatchException - * if the matrices dimensions do not match. - * @throws SingularMatrixException - * if the decomposed matrix is singular. - */ - FieldVector<T> solve(final FieldVector<T> b); - - /** Solve the linear equation A × X = B for matrices A. - * <p>The A matrix is implicit, it is provided by the underlying - * decomposition algorithm.</p> - * @param b right-hand side of the equation A × X = B - * @return a matrix X that minimizes the two norm of A × X - B - * @throws org.apache.commons.math3.exception.DimensionMismatchException - * if the matrices dimensions do not match. - * @throws SingularMatrixException - * if the decomposed matrix is singular. - */ - FieldMatrix<T> solve(final FieldMatrix<T> b); - - /** - * Check if the decomposed matrix is non-singular. - * @return true if the decomposed matrix is non-singular - */ - boolean isNonSingular(); - - /** Get the inverse (or pseudo-inverse) of the decomposed matrix. - * @return inverse matrix - * @throws SingularMatrixException - * if the decomposed matrix is singular. - */ - FieldMatrix<T> getInverse(); -}
