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 imagingbook.common.math.Matrix; 012import org.apache.commons.math3.linear.MatrixUtils; 013import org.apache.commons.math3.linear.RealMatrix; 014import org.apache.commons.math3.linear.RealVector; 015 016import static imagingbook.common.math.Matrix.isSquare; 017import static imagingbook.common.math.Matrix.sameSize; 018import static imagingbook.common.math.eigen.eispack.QZHES.qzhes; 019import static imagingbook.common.math.eigen.eispack.QZIT.qzit; 020import static imagingbook.common.math.eigen.eispack.QZVAL.qzval; 021import static imagingbook.common.math.eigen.eispack.QZVEC.qzvec; 022 023/** 024 * <p> 025 * Solves the generalized eigenproblem of the form A x = λ B x, for square matrices A, B. 026 * </p> 027 * <p> 028 * This implementation was ported from original EISPACK code [1] using a finite state machine concept [2] to untangle 029 * Fortran's GOTO statements (which are not available in Java), with some inspirations from this <a href= 030 * "https://github.com/accord-net/framework/blob/development/Sources/Accord.Math/Decompositions/GeneralizedEigenvalueDecomposition.cs"> 031 * C# implementation</a>. See 032 * <a href="https://mipav.cit.nih.gov/">mipav.cit.nih.gov</a> (file 033 * {@code gov.nih.mipav.model.structures.jama.GeneralizedEigenvalue.java}) for another Java implementation based on 034 * EISPACK. 035 * </p> 036 * <p> 037 * Note: Results have limited accuracy. For some reason the first eigenvalue/-vector (k=0) is reasonably accurate, but 038 * the remaining ones are not. If matrices A, B are symmetric and B is positive definite, better use 039 * {@link GeneralizedSymmetricEigenDecomposition} instead, which is more accurate. 040 * </p> 041 * <p> 042 * [1] http://www.netlib.no/netlib/eispack/ <br> [2] D. E. Knuth, "Structured Programming with Goto Statements", 043 * Computing Surveys, Vol. 6, No. 4 (1974). 044 * 045 * @author WB 046 * @version 2022/12/22 047 * @see GeneralizedSymmetricEigenDecomposition 048 */ 049public class GeneralizedEigenDecomposition implements RealEigenDecomposition { 050 // TODO: Check complex-valued eigenvectors. 051 052 public static boolean VERBOSE = false; 053 054 private final int n; 055 private final double[] alphaR; 056 private final double[] alphaI; 057 private final double[] beta; 058 private final double[][] Z; 059 060 /** 061 * Constructor, solves the generalized eigenproblem A x = λ B x, for square matrices A, B. 062 * 063 * @param A square matrix A 064 * @param B square matrix B 065 */ 066 public GeneralizedEigenDecomposition(RealMatrix A, RealMatrix B) { 067 this(A.getData(), B.getData()); 068 } 069 070 /** 071 * Constructor, solves the generalized eigenproblem A x = λ B x, for square matrices A, B. 072 * 073 * @param A square matrix A 074 * @param B square matrix B 075 */ 076 public GeneralizedEigenDecomposition(double[][] A, double[][] B) { 077 if (!isSquare(A) || !isSquare(B) || !sameSize(A, B)) { 078 throw new IllegalArgumentException("matrices A, B must be square and of same size"); 079 } 080 this.n = A.length; 081 this.alphaR = new double[n]; 082 this.alphaI = new double[n]; 083 this.beta = new double[n]; 084 this.Z = new double[n][n]; 085 086 // a, b are modified by EISPACK routines: 087 double[][] a = Matrix.duplicate(A); 088 double[][] b = Matrix.duplicate(B); 089 boolean matz = true; 090 091 qzhes(a, b, matz, Z); 092 int ierr = qzit(a, b, 0, matz, Z); 093 if (ierr >= 0) { 094 throw new RuntimeException("limit of 30*n iterations exhausted for eigenvalue " + ierr); 095 } 096 qzval(a, b, alphaR, alphaI, beta, matz, Z); 097 qzvec(a, b, alphaR, alphaI, beta, Z); 098 } 099 100 /** 101 * Returns a vector with the real parts of the eigenvalues. 102 * 103 * @return the real parts of the eigenvalues 104 */ 105 @Override 106 public double[] getRealEigenvalues() { 107 double[] eval = new double[n]; 108 for (int i = 0; i < n; i++) { 109 eval[i] = alphaR[i] / beta[i]; 110 } 111 return eval; 112 } 113 114 @Override 115 public double getRealEigenvalue(int k) { 116 return alphaR[k] / beta[k]; 117 } 118 119 /** 120 * Returns a vector with the imaginary parts of the eigenvalues. 121 * 122 * @return the imaginary parts of the eigenvalues 123 */ 124 public double[] getImagEigenvalues() { 125 double[] eval = new double[n]; 126 for (int i = 0; i < n; i++) 127 eval[i] = alphaI[i] / beta[i]; 128 return eval; 129 } 130 131 /** 132 * Return the matrix of eigenvectors, which are its column vectors. 133 * 134 * @return the matrix of eigenvectors 135 */ 136 @Override 137 public RealMatrix getV() { // TODO: not implemented/checked for complex eigenvalues yet! 138 return MatrixUtils.createRealMatrix(Z); 139 } 140 141 /** 142 * Returns the specified eigenvector. 143 * 144 * @param k index 145 * @return the kth eigenvector 146 */ 147 @Override 148 public RealVector getEigenvector(int k) { // TODO: not implemented for complex eigenvalues yet! 149 return MatrixUtils.createRealVector(Matrix.getColumn(Z, k)); 150 } 151 152 /** 153 * Returns whether the calculated eigenvalues are complex or real. The method performs a zero check on each element 154 * of the {@link #getImagEigenvalues()} array and returns {@code true} if any element is not equal to zero. 155 * 156 * @return {@code true} if any of the eigenvalues is complex, {@code false} otherwise 157 */ 158 @Override 159 public boolean hasComplexEigenvalues() { 160 for (int i = 0; i < n; i++) { 161 if (alphaI[i] != 0.0) { 162 return true; 163 } 164 } 165 return false; 166 } 167 168 /** 169 * Returns the block diagonal eigenvalue matrix D. 170 * 171 * @return the block diagonal eigenvalue matrix 172 */ 173 @Override 174 public RealMatrix getD() { // TODO: check complex case! 175 double[][] x = new double[n][n]; 176 177 for (int i = 0; i < n; i++) { 178 for (int j = 0; j < n; j++) 179 x[i][j] = 0.0; 180 181 x[i][i] = alphaR[i] / beta[i]; 182 if (alphaI[i] > 0) 183 x[i][i + 1] = alphaI[i] / beta[i]; 184 else if (alphaI[i] < 0) 185 x[i][i - 1] = alphaI[i] / beta[i]; 186 } 187 188 return MatrixUtils.createRealMatrix(x); 189 } 190 191 // ---------------------------------------------------------- 192 193 /** 194 * This method checks if any of the generated betas is zero. It does not says that the problem is singular, but only 195 * that one of the matrices A, B is singular. 196 * 197 * @return true if A or B is singular 198 */ 199 public boolean isSingular() { 200 for (int i = 0; i < n; i++) { 201 if (beta[i] == 0) 202 return true; 203 } 204 return false; 205 } 206 207 /** 208 * Returns true if the eigenvalue problem is degenerate (i.e., ill-posed). 209 * 210 * @return true if degenerate 211 */ 212 public boolean IsDegenerate() { 213 for (int i = 0; i < n; i++) { 214 if (beta[i] == 0 && alphaR[i] == 0) 215 return true; 216 } 217 return false; 218 } 219 220}