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 ******************************************************************************/ 009 010package imagingbook.common.math.eigen; 011 012import imagingbook.common.math.Matrix; 013import org.apache.commons.math3.linear.MatrixUtils; 014import org.apache.commons.math3.linear.RealMatrix; 015import org.apache.commons.math3.linear.RealVector; 016 017import static imagingbook.common.math.Arithmetic.sqr; 018 019/** 020 * <p> 021 * Implements an efficient, closed form algorithm for calculating the real eigenvalues (λ) and eigenvectors (x) 022 * of a 2x2 matrix of the form 023 * </p> 024 * <pre> 025 * | a b | 026 * | c d | 027 * </pre> 028 * <p> 029 * There are typically (but not always) two pairs of real-valued solutions ⟨λ<sub>0</sub>, 030 * x<sub>0</sub>⟩, ⟨λ<sub>1</sub>, x<sub>1</sub>⟩ such that A·x<sub>k</sub> = 031 * λ<sub>k</sub>·x<sub>k</sub>. The resulting eigensystems are ordered such that λ<sub>0</sub> ≥ 032 * λ<sub>1</sub>. Eigenvectors are not normalized, i.e., no unit vectors (any scalar multiple of an Eigenvector 033 * is an Eigenvector too). Non-real eigenvalues are not handled. Clients should call method 034 * {@link #hasComplexEigenvalues()} to check if the eigenvalue calculation was successful. 035 * </p> 036 * <p> 037 * This implementation is inspired by Ch. 5 ("Consider the Lowly 2x2 Matrix") of [1]. Note that Blinn uses the notation 038 * x·A = λ·x for the matrix-vector product (as common in computer graphics), while this 039 * implementation uses A·x = λ·x. Thus x is treated as a column vector and matrix A is transposed 040 * (elements b/c are exchanged). See Appendix Sec. B.5 of [2] for more details. 041 * </p> 042 * <p> 043 * This implementation is considerably faster (ca. factor 5) than the general solution (e.g., 044 * {@link EigenDecompositionJama}) for 2x2 matrices. 045 * </p> 046 * <p> 047 * [1] Blinn, Jim: "Jim Blinn's Corner: Notation, Notation, Notation", Morgan Kaufmann (2002). <br> [2] W. Burger, M.J. 048 * Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer (2022). 049 * </p> 050 * 051 * @author WB 052 * @version 2022/02/18 053 * @see EigenDecompositionJama 054 * @see EigenDecompositionApache 055 */ 056public class Eigensolver2x2 implements RealEigenDecomposition { // to check: http://www.akiti.ca/Eig2Solv.html 057 058 private final boolean isReal; 059 private final double[] eVals = {Double.NaN, Double.NaN}; 060 private final double[][] eVecs = new double[2][]; 061 062 /** 063 * Constructor, takes a 2x2 matrix. 064 * @param A a 2x2 matrix 065 */ 066 public Eigensolver2x2(double[][] A) { 067 this(A[0][0], A[0][1], A[1][0], A[1][1]); 068 if (Matrix.getNumberOfRows(A) != 2 || Matrix.getNumberOfColumns(A) != 2) { 069 throw new IllegalArgumentException("matrix not of size 2x2"); 070 } 071 } 072 073 /** 074 * Constructor, takes the individual elements of a 2x2 matrix A: 075 * <pre> 076 * | a b | 077 * | c d | </pre> 078 * @param a matrix element A[0,0] 079 * @param b matrix element A[0,1] 080 * @param c matrix element A[1,0] 081 * @param d matrix element A[1,1] 082 */ 083 public Eigensolver2x2(double a, double b, double c, double d) { 084 isReal = solve(a, b, c, d); 085 } 086 087 private boolean solve(final double a, final double b, final double c, final double d) { 088 final double R = (a + d) / 2; 089 final double S = (a - d) / 2; 090 final double rho = sqr(S) + b * c; 091 092 if (rho < 0) { 093 return false; // no real-valued eigenvalues 094 } 095 096 final double T = Math.sqrt(rho); 097 final double lambda0 = R + T; // eigenvalue 0 098 final double lambda1 = R - T; // eigenvalue 1 099 final double[] x0, x1; // eigenvectors 0, 1 100 101 if (a - d > 0) { 102// System.out.println("Case 1"); 103 x0 = new double[] {S + T, c}; 104 x1 = new double[] {b, -(S + T)}; 105 } 106 else if (a - d < 0) { 107// System.out.println("Case 2"); 108 x0 = new double[] {b, -S + T}; 109 x1 = new double[] {S - T, c}; 110 } 111 else { // (A - D) == 0 112// System.out.println("Case 3"); 113 final double bA = Math.abs(b); 114 final double cA = Math.abs(c); 115 final double bcR = Math.sqrt(b * c); 116 if (bA < cA) { // |b| < |c| 117// System.out.println("Case 3a"); 118 x0 = new double[] { bcR, c}; 119 x1 = new double[] {-bcR, c}; 120 } 121 else if (bA > cA) { // |b| > |c| 122// System.out.println("Case 3b"); 123 x0 = new double[] {b, bcR}; 124 x1 = new double[] {b, -bcR}; 125 } 126 else { // |B| == |C| and B,C must have the same sign 127// System.out.println("Case 3c"); 128 if (cA > 0) { // 129 x0 = new double[] { c, c}; 130 x1 = new double[] {-c, c}; 131 } 132 else { // B = C = 0; any vector is an eigenvector (we don't return trivial zero vectors) 133 x0 = new double[] { 0, 1}; // pick 2 arbitrary, orthogonal vectors 134 x1 = new double[] { 1, 0}; 135 } 136 } 137 } 138 139 eVals[0] = lambda0; 140 eVals[1] = lambda1; 141 eVecs[0] = x0; 142 eVecs[1] = x1; 143 144 // lambda0 >= lambda1, no need to sort by magnitude 145 146// if (Math.abs(lambda0) >= Math.abs(lambda1)) { // order eigenvalues by magnitude 147// eVals[0] = lambda0; 148// eVals[1] = lambda1; 149// eVecs[0] = x0; 150// eVecs[1] = x1; 151// } 152// else { 153// eVals[0] = lambda1; 154// eVals[1] = lambda0; 155// eVecs[0] = x1; 156// eVecs[1] = x0; 157// } 158 return true; // real eigenvalues 159 } 160 161 162// public boolean isReal() { 163// return isReal; 164// } 165 166 @Override 167 public boolean hasComplexEigenvalues() { 168 return !isReal; 169 } 170 171 @Override 172 public double[] getRealEigenvalues() { 173 return eVals; 174 } 175 176 @Override 177 public double getRealEigenvalue(int k) { 178 return eVals[k]; 179 } 180 181 @Override 182 public RealMatrix getV() { 183 return MatrixUtils.createRealMatrix(Matrix.transpose(eVecs)); 184 } 185 186 @Override 187 public RealVector getEigenvector(int k) { 188 return MatrixUtils.createRealVector(eVecs[k]); 189 } 190 191 @Override 192 public String toString() { 193 if (this.isReal) { 194 return String.format("<%.4f, %.4f, %s, %s>", 195 eVals[0], eVals[1], Matrix.toString(eVecs[0]), Matrix.toString(eVecs[1])); 196 } 197 else { 198 return "<not real>"; 199 } 200 } 201 202 203}