001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 * image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause).
006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved.
007 * Visit https://imagingbook.com for additional details.
008 ******************************************************************************/
009package imagingbook.common.math.eigen;
010
011import org.apache.commons.math3.linear.CholeskyDecomposition;
012import org.apache.commons.math3.linear.DecompositionSolver;
013import org.apache.commons.math3.linear.EigenDecomposition;
014import org.apache.commons.math3.linear.LUDecomposition;
015import org.apache.commons.math3.linear.MatrixUtils;
016import org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException;
017import org.apache.commons.math3.linear.RealMatrix;
018import org.apache.commons.math3.linear.RealVector;
019
020import static org.apache.commons.math3.linear.CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD;
021import static org.apache.commons.math3.linear.CholeskyDecomposition.DEFAULT_RELATIVE_SYMMETRY_THRESHOLD;
022
023/**
024 * <p>
025 * Solves the generalized symmetric eigenproblem of the form A x = &lambda; B x, where matrices A, B are symmetric and B
026 * is positive definite (see Sec. 11.0.5. of [1]). See Appendix Sec. B.5.2 of [2] for more details. The methods defined
027 * by this class are analogous to the conventional eigendecomposition (see {@link EigenDecomposition}).
028 * </p>
029 * <p>
030 * [1] Press, Teukolsky, Vetterling, Flannery: "Numerical Recipes". Cambridge University Press, 3rd ed. (2007). <br> [2]
031 * W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
032 * (2022).
033 * </p>
034 *
035 * @author WB
036 * @version 2022/06/11
037 * @see GeneralizedEigenDecomposition
038 * @see EigenDecompositionJama
039 */
040public class GeneralizedSymmetricEigenDecomposition implements RealEigenDecomposition {
041
042        private final EigenDecomposition eigendecompY;
043        private final DecompositionSolver solverLT;
044
045        /**
046         * Constructor, using specific parameters. Solves the generalized symmetric eigenproblem of the form A x = &lambda;
047         * B x, where matrices A, B are symmetric and B is positive definite. An exception is thrown if A is not symmetric
048         * and the Cholesky decomposition throws an exception if B is either not symmetric, not positive definite or
049         * singular.
050         *
051         * @param A real symmetric matrix
052         * @param B real symmetric and positive definite matrix
053         * @param rsth relative symmetry threshold
054         * @param apth absolute positivity threshold
055         */
056        public GeneralizedSymmetricEigenDecomposition(RealMatrix A, RealMatrix B, double rsth, double apth) {
057                if (!MatrixUtils.isSymmetric(A, rsth)) {
058                        throw new IllegalArgumentException("matrix A must be symmetric");
059                }
060
061                if (!MatrixUtils.isSymmetric(B, rsth)) {
062                        throw new IllegalArgumentException("matrix B must be symmetric");
063                }
064
065                CholeskyDecomposition cd = null;
066                try {
067                        cd = new CholeskyDecomposition(B, rsth, apth);
068                } catch (NonPositiveDefiniteMatrixException e) {
069                        throw new IllegalArgumentException("matrix B must be positive definite");
070                }
071
072                RealMatrix L = cd.getL();
073
074                // find Q, such that Q * LT = A or equivalently L * QT = AT = A (since A is symmetric)
075                DecompositionSolver sL = new LUDecomposition(L).getSolver();
076                RealMatrix Q = sL.solve(A);
077
078                // find Y, such that L * Y = QT
079                RealMatrix Y = sL.solve(Q.transpose());
080
081                // Y has the same eigenvalues as the original system and eigenvectors v_k
082                this.eigendecompY = new EigenDecomposition(Y);    // use EigenDecompositionJama instead?
083
084                // the eigenvectors x_k of the original system are related
085                // to the eigenvectors v_k of Y as x_k = (LT)^(-1) * v_k = (L^(-1))^T * v_k
086                // or found by solving L^T * x_k = y_k, using the following solver:
087                this.solverLT = new LUDecomposition(cd.getLT()).getSolver();
088        }
089
090        /**
091         * Constructor, using default parameters. Solves the generalized symmetric eigenproblem of the form A x = &lambda; B
092         * x, where matrices A, B * are symmetric and B is positive definite
093         *
094         * @param A real symmetric matrix
095         * @param B real symmetric and positive definite matrix
096         */
097        public GeneralizedSymmetricEigenDecomposition(RealMatrix A, RealMatrix B) {
098                this(A, B,
099                                DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
100                                DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
101        }
102
103        // ---------------------------------------------------------------------
104
105        @Override
106        public double[] getRealEigenvalues() {
107                return eigendecompY.getRealEigenvalues();
108        }
109
110        @Override
111        public double getRealEigenvalue(int k) {
112                return eigendecompY.getRealEigenvalue(k);
113        }
114
115        public double[] getImagEigenvalues() {
116                return eigendecompY.getImagEigenvalues();
117        }
118
119        @Override
120        public boolean hasComplexEigenvalues() {
121                return eigendecompY.hasComplexEigenvalues();
122        }
123
124        @Override
125        public RealMatrix getD() {
126                return eigendecompY.getD();
127        }
128
129        @Override
130        public RealVector getEigenvector(int k) {
131//              return LiT.operate(ed.getEigenvector(k));
132                // solve LT * x_k = v_k
133                RealVector vk = eigendecompY.getEigenvector(k);
134                return solverLT.solve(vk);
135        }
136
137        @Override
138        public RealMatrix getV() {
139//              return LiT.multiply(ed.getV());
140                // solve LT * X = V
141                return solverLT.solve(eigendecompY.getV());
142        }
143
144}
145