The eigen decomposition of matrix A is a set of two matrices: - * V and D such that A = V D VT. A, V and D are all m × m - * matrices.
- *As of 2.0, this class supports only symmetric matrices, - * and hence computes only real realEigenvalues. This implies the D matrix returned by - * {@link #getD()} is always diagonal and the imaginary values returned {@link - * #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always null.
- *When called with a {@link RealMatrix} argument, this implementation only uses - * the upper part of the matrix, the part below the diagonal is not accessed at all.
- *Eigenvalues are computed as soon as the matrix is decomposed, but eigenvectors - * are computed only when required, i.e. only when one of the {@link #getEigenvector(int)}, - * {@link #getV()}, {@link #getVT()}, {@link #getSolver()} methods is called.
- *This implementation is based on Inderjit Singh Dhillon thesis - * A - * New O(n2) Algorithm for the Symmetric Tridiagonal Eigenvalue/Eigenvector - * Problem, on Beresford N. Parlett and Osni A. Marques paper An Implementation of the - * dqds Algorithm (Positive Case) and on the corresponding LAPACK routines (DLARRE, - * DLASQ2, DLAZQ3, DLAZQ4, DLASQ5 and DLASQ6).
- *The authors of the original fortran version are: - *
+ * The eigen decomposition of matrix A is a set of two matrices: V and D such + * that A = V D VT. A, V and D are all m × m matrices. + *
+ *+ * As of 2.0, this class supports only symmetric matrices, and + * hence computes only real realEigenvalues. This implies the D matrix returned + * by {@link #getD()} is always diagonal and the imaginary values returned + * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always + * null. + *
+ *+ * When called with a {@link RealMatrix} argument, this implementation only uses + * the upper part of the matrix, the part below the diagonal is not accessed at + * all. + *
+ *+ * 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 *
* @version $Revision$ $Date$ * @since 2.0 */ public class EigenDecompositionImpl implements EigenDecomposition { - /** Tolerance. */ - private static final double TOLERANCE = 100 * MathUtils.EPSILON; - - /** Squared tolerance. */ - private static final double TOLERANCE_2 = TOLERANCE * TOLERANCE; - - /** Split tolerance. */ - private double splitTolerance; + /** Maximum number of iterations accepted in the implicit QL transformation */ + private byte maxIter = 30; /** Main diagonal of the tridiagonal matrix. */ private double[] main; @@ -75,66 +60,12 @@ public class EigenDecompositionImpl implements EigenDecomposition { /** Secondary diagonal of the tridiagonal matrix. */ private double[] secondary; - /** Squared secondary diagonal of the tridiagonal matrix. */ - private double[] squaredSecondary; - - /** Transformer to tridiagonal (may be null if matrix is already tridiagonal). */ + /** + * Transformer to tridiagonal (may be null if matrix is already + * tridiagonal). + */ private TriDiagonalTransformer transformer; - /** Lower bound of spectra. */ - private double lowerSpectra; - - /** Upper bound of spectra. */ - private double upperSpectra; - - /** Minimum pivot in the Sturm sequence. */ - private double minPivot; - - /** Current shift. */ - private double sigma; - - /** Low part of the current shift. */ - private double sigmaLow; - - /** Shift increment to apply. */ - private double tau; - - /** Work array for all decomposition algorithms. */ - private double[] work; - - /** Shift within qd array for ping-pong implementation. */ - private int pingPong; - - /** Max value of diagonal elements in current segment. */ - private double qMax; - - /** Min value of off-diagonal elements in current segment. */ - private double eMin; - - /** Type of the last dqds shift. */ - private int tType; - - /** Minimal value on current state of the diagonal. */ - private double dMin; - - /** Minimal value on current state of the diagonal, excluding last element. */ - private double dMin1; - - /** Minimal value on current state of the diagonal, excluding last two elements. */ - private double dMin2; - - /** Last value on current state of the diagonal. */ - private double dN; - - /** Last but one value on current state of the diagonal. */ - private double dN1; - - /** Last but two on current state of the diagonal. */ - private double dN2; - - /** Shift ratio with respect to dMin used when tType == 6. */ - private double g; - /** Real part of the realEigenvalues. */ private double[] realEigenvalues; @@ -156,70 +87,65 @@ public class EigenDecompositionImpl implements EigenDecomposition { /** * Calculates the eigen decomposition of the given symmetric matrix. * @param matrix The symmetric matrix to decompose. - * @param splitTolerance tolerance on the off-diagonal elements relative to the - * geometric mean to split the tridiagonal matrix (a suggested value is - * {@link MathUtils#SAFE_MIN}) - * @exception InvalidMatrixException (wrapping a {@link - * org.apache.commons.math.ConvergenceException} if algorithm fails to converge + * @param splitTolerance dummy parameter, present for backward compatibility only. + * @exception InvalidMatrixException (wrapping a + * {@link org.apache.commons.math.ConvergenceException} if algorithm + * fails to converge */ - public EigenDecompositionImpl(final RealMatrix matrix, - final double splitTolerance) - throws InvalidMatrixException { + public EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance) + throws InvalidMatrixException { if (isSymmetric(matrix)) { - this.splitTolerance = splitTolerance; transformToTridiagonal(matrix); - decompose(); + findEigenVectors(transformer.getQ().getData()); } else { - // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are NOT supported + // as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are + // NOT supported // see issue https://issues.apache.org/jira/browse/MATH-235 - throw new InvalidMatrixException("eigen decomposition of assymetric matrices not supported yet"); + throw new InvalidMatrixException( + "eigen decomposition of assymetric matrices not supported yet"); } } /** - * Calculates the eigen decomposition of the given tridiagonal symmetric matrix. - * @param main the main diagonal of the matrix (will be copied) - * @param secondary the secondary diagonal of the matrix (will be copied) - * @param splitTolerance tolerance on the off-diagonal elements relative to the - * geometric mean to split the tridiagonal matrix (a suggested value is - * {@link MathUtils#SAFE_MIN}) - * @exception InvalidMatrixException (wrapping a {@link - * org.apache.commons.math.ConvergenceException} if algorithm fails to converge + * 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 triadiagonal form + * @param secondary Secondary of the tridiagonal form + * @param splitTolerance dummy parameter, present for backward compatibility only. + * @exception InvalidMatrixException (wrapping a + * {@link org.apache.commons.math.ConvergenceException} if algorithm + * fails to converge */ - public EigenDecompositionImpl(final double[] main, double[] secondary, + public EigenDecompositionImpl(final double[] main,final double[] secondary, final double splitTolerance) - throws InvalidMatrixException { - + throws InvalidMatrixException { this.main = main.clone(); this.secondary = secondary.clone(); transformer = null; - - // pre-compute some elements - squaredSecondary = new double[secondary.length]; - for (int i = 0; i < squaredSecondary.length; ++i) { - final double s = secondary[i]; - squaredSecondary[i] = s * s; + final int size=main.length; + double[][] z = new double[size][size]; + for (int i=0;iThis method only find exact linear solutions, i.e. solutions for - * which ||A × X - B|| is exactly 0.
- * @param b right-hand side of the equation A × X = B + /** + * Solve the linear equation A × X = B for symmetric matrices A. + *+ * This method only find exact linear solutions, i.e. solutions for + * which ||A × X - B|| is exactly 0. + *
+ * @param b + * right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B - * @exception IllegalArgumentException if matrices dimensions don't match - * @exception InvalidMatrixException if decomposed matrix is singular + * @exception IllegalArgumentException + * if matrices dimensions don't match + * @exception InvalidMatrixException + * if decomposed matrix is singular */ public double[] solve(final double[] b) - throws IllegalArgumentException, InvalidMatrixException { + throws IllegalArgumentException, InvalidMatrixException { if (!isNonSingular()) { throw new SingularMatrixException(); @@ -420,16 +309,22 @@ public class EigenDecompositionImpl implements EigenDecomposition { } - /** Solve the linear equation A × X = B for symmetric matrices A. - *This method only find exact linear solutions, i.e. solutions for - * which ||A × X - B|| is exactly 0.
- * @param b right-hand side of the equation A × X = B + /** + * Solve the linear equation A × X = B for symmetric matrices A. + *+ * This method only find exact linear solutions, i.e. solutions for + * which ||A × X - B|| is exactly 0. + *
+ * @param b + * right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B - * @exception IllegalArgumentException if matrices dimensions don't match - * @exception InvalidMatrixException if decomposed matrix is singular + * @exception IllegalArgumentException + * if matrices dimensions don't match + * @exception InvalidMatrixException + * if decomposed matrix is singular */ public RealVector solve(final RealVector b) - throws IllegalArgumentException, InvalidMatrixException { + throws IllegalArgumentException, InvalidMatrixException { if (!isNonSingular()) { throw new SingularMatrixException(); @@ -438,8 +333,8 @@ public class EigenDecompositionImpl implements EigenDecomposition { final int m = realEigenvalues.length; if (b.getDimension() != m) { throw MathRuntimeException.createIllegalArgumentException( - "vector length mismatch: got {0} but expected {1}", - b.getDimension(), m); + "vector length mismatch: got {0} but expected {1}", b + .getDimension(), m); } final double[] bp = new double[m]; @@ -456,16 +351,22 @@ public class EigenDecompositionImpl implements EigenDecomposition { } - /** Solve the linear equation A × X = B for symmetric matrices A. - *This method only find exact linear solutions, i.e. solutions for - * which ||A × X - B|| is exactly 0.
- * @param b right-hand side of the equation A × X = B + /** + * Solve the linear equation A × X = B for symmetric matrices A. + *+ * This method only find exact linear solutions, i.e. solutions for + * which ||A × X - B|| is exactly 0. + *
+ * @param b + * right-hand side of the equation A × X = B * @return a matrix X that minimizes the two norm of A × X - B - * @exception IllegalArgumentException if matrices dimensions don't match - * @exception InvalidMatrixException if decomposed matrix is singular + * @exception IllegalArgumentException + * if matrices dimensions don't match + * @exception InvalidMatrixException + * if decomposed matrix is singular */ public RealMatrix solve(final RealMatrix b) - throws IllegalArgumentException, InvalidMatrixException { + throws IllegalArgumentException, InvalidMatrixException { if (!isNonSingular()) { throw new SingularMatrixException(); @@ -473,9 +374,11 @@ public class EigenDecompositionImpl implements EigenDecomposition { final int m = realEigenvalues.length; if (b.getRowDimension() != m) { - throw MathRuntimeException.createIllegalArgumentException( - "dimensions mismatch: got {0}x{1} but expected {2}x{3}", - b.getRowDimension(), b.getColumnDimension(), m, "n"); + throw MathRuntimeException + .createIllegalArgumentException( + "dimensions mismatch: got {0}x{1} but expected {2}x{3}", + b.getRowDimension(), b.getColumnDimension(), m, + "n"); } final int nColB = b.getColumnDimension(); @@ -512,12 +415,13 @@ public class EigenDecompositionImpl implements EigenDecomposition { return true; } - /** Get the inverse of the decomposed matrix. + /** + * Get the inverse of the decomposed matrix. * @return inverse matrix - * @throws InvalidMatrixException if decomposed matrix is singular + * @throws InvalidMatrixException + * if decomposed matrix is singular */ - public RealMatrix getInverse() - throws InvalidMatrixException { + public RealMatrix getInverse() throws InvalidMatrixException { if (!isNonSingular()) { throw new SingularMatrixException(); @@ -545,1382 +449,146 @@ public class EigenDecompositionImpl implements EigenDecomposition { /** * Transform matrix to tridiagonal. - * @param matrix matrix to transform + * @param matrix + * matrix to transform */ private void transformToTridiagonal(final RealMatrix matrix) { // transform the matrix to tridiagonal transformer = new TriDiagonalTransformer(matrix); - main = transformer.getMainDiagonalRef(); + main = transformer.getMainDiagonalRef(); secondary = transformer.getSecondaryDiagonalRef(); - // pre-compute some elements - squaredSecondary = new double[secondary.length]; - for (int i = 0; i < squaredSecondary.length; ++i) { - final double s = secondary[i]; - squaredSecondary[i] = s * s; - } - } /** - * Compute the Gershgorin circles for all rows. + * Find eigenvalues and eigenvectors (Dubrulle et al., 1971) + * @param householderMatrix Householder matrix of the transformation + * to tri-diagonal form. */ - private void computeGershgorinCircles() { - - final int m = main.length; - final int lowerStart = 4 * m; - final int upperStart = 5 * m; - lowerSpectra = Double.POSITIVE_INFINITY; - upperSpectra = Double.NEGATIVE_INFINITY; - double eMax = 0; - - double eCurrent = 0; - for (int i = 0; i < m - 1; ++i) { - - final double dCurrent = main[i]; - final double ePrevious = eCurrent; - eCurrent = Math.abs(secondary[i]); - eMax = Math.max(eMax, eCurrent); - final double radius = ePrevious + eCurrent; - - final double lower = dCurrent - radius; - work[lowerStart + i] = lower; - lowerSpectra = Math.min(lowerSpectra, lower); - - final double upper = dCurrent + radius; - work[upperStart + i] = upper; - upperSpectra = Math.max(upperSpectra, upper); + private void findEigenVectors(double[][] householderMatrix) { + double[][]z = householderMatrix.clone(); + final int n = main.length; + realEigenvalues = new double[n]; + imagEigenvalues = new double[n]; + double[] e = new double[n]; + for (int i = 0; i < n - 1; i++) { + realEigenvalues[i] = main[i]; + e[i] = secondary[i]; } - - final double dCurrent = main[m - 1]; - final double lower = dCurrent - eCurrent; - work[lowerStart + m - 1] = lower; - lowerSpectra = Math.min(lowerSpectra, lower); - final double upper = dCurrent + eCurrent; - work[upperStart + m - 1] = upper; - upperSpectra = Math.max(upperSpectra, upper); - minPivot = MathUtils.SAFE_MIN * Math.max(1.0, eMax * eMax); - - } - - /** - * Find the realEigenvalues. - * @exception InvalidMatrixException if a block cannot be diagonalized - */ - private void findEigenvalues() - throws InvalidMatrixException { - - // compute splitting points - ListIn low dimensions, we simply solve the characteristic polynomial.
- * @param index index of the first row of the block - */ - private void process1RowBlock(final int index) { - realEigenvalues[index] = main[index]; - } - - /** - * Find realEigenvalues in a block with 2 rows. - *In low dimensions, we simply solve the characteristic polynomial.
- * @param index index of the first row of the block - * @exception InvalidMatrixException if characteristic polynomial cannot be solved - */ - private void process2RowsBlock(final int index) - throws InvalidMatrixException { - - // the characteristic polynomial is - // X^2 - (q0 + q1) X + q0 q1 - e1^2 - final double q0 = main[index]; - final double q1 = main[index + 1]; - final double e12 = squaredSecondary[index]; - - final double s = q0 + q1; - final double p = q0 * q1 - e12; - final double delta = s * s - 4 * p; - if (delta < 0) { - throw new InvalidMatrixException("cannot solve degree {0} equation", 2); - } - - final double largestRoot = 0.5 * (s + Math.sqrt(delta)); - realEigenvalues[index] = largestRoot; - realEigenvalues[index + 1] = p / largestRoot; - - } - - /** - * Find realEigenvalues in a block with 3 rows. - *In low dimensions, we simply solve the characteristic polynomial.
- * @param index index of the first row of the block - * @exception InvalidMatrixException if diagonal elements are not positive - */ - private void process3RowsBlock(final int index) - throws InvalidMatrixException { - - // the characteristic polynomial is - // X^3 - (q0 + q1 + q2) X^2 + (q0 q1 + q0 q2 + q1 q2 - e1^2 - e2^2) X + q0 e2^2 + q2 e1^2 - q0 q1 q2 - final double q0 = main[index]; - final double q1 = main[index + 1]; - final double q2 = main[index + 2]; - final double e12 = squaredSecondary[index]; - final double q1q2Me22 = q1 * q2 - squaredSecondary[index + 1]; - - // compute coefficients of the cubic equation as: x^3 + b x^2 + c x + d = 0 - final double b = -(q0 + q1 + q2); - final double c = q0 * q1 + q0 * q2 + q1q2Me22 - e12; - final double d = q2 * e12 - q0 * q1q2Me22; - - // solve cubic equation - final double b2 = b * b; - final double beta = b / 3; - final double q = (3 * c - b2) / 9; - final double r = ((9 * c - 2 * b2) * b - 27 * d) / 54; - final double delta = q * q * q + r * r; - double z0; - double z1; - double z2; - if (delta > 0) { - // there are two complex solutions, we cannot handle this - throw new InvalidMatrixException("cannot solve degree {0} equation", 3); - } else if (delta < 0) { - // three different real roots - final double sqrtMq = Math.sqrt(-q); - final double theta = Math.acos(r / (-q * sqrtMq)); - final double alpha = 2 * sqrtMq; - z0 = alpha * Math.cos(theta / 3) - beta; - z1 = alpha * Math.cos((theta + 2 * Math.PI) / 3) - beta; - z2 = alpha * Math.cos((theta + 4 * Math.PI) / 3) - beta; - } else { - // three real roots, two of which are equal - final double cbrtR = Math.cbrt(r); - z0 = 2 * cbrtR - beta; - z1 = -cbrtR - beta; - z2 = z1; - } - - // sort the eigenvalues - if (z0 < z1) { - final double t = z0; - z0 = z1; - z1 = t; - } - if (z1 < z2) { - final double t = z1; - z1 = z2; - z2 = t; - } - if (z0 < z1) { - final double t = z0; - z0 = z1; - z1 = t; - } - - realEigenvalues[index] = z0; - realEigenvalues[index + 1] = z1; - realEigenvalues[index + 2] = z2; - - } - - /** - * Find realEigenvalues using dqd/dqds algorithms. - *This implementation is based on Beresford N. Parlett - * and Osni A. Marques paper An - * Implementation of the dqds Algorithm (Positive Case) and on the - * corresponding LAPACK routine DLASQ2.
- * @param n number of rows of the block - * @exception InvalidMatrixException if block cannot be diagonalized - * after 30 * n iterations - */ - private void processGeneralBlock(final int n) - throws InvalidMatrixException { - - // check decomposed matrix data range - double sumOffDiag = 0; - for (int i = 0; i < n - 1; ++i) { - final int fourI = 4 * i; - final double ei = work[fourI + 2]; - sumOffDiag += ei; - } - - if (sumOffDiag == 0) { - // matrix is already diagonal - return; - } - - // initial checks for splits (see Parlett & Marques section 3.3) - flipEveryOtherIfWarranted(n); - - // two iterations with Li's test for initial splits - initialSplits(n); - - // initialize parameters used by goodStep - tType = 0; - dMin1 = 0; - dMin2 = 0; - dN = 0; - dN1 = 0; - dN2 = 0; - tau = 0; - - // process split segments - int i0 = 0; - int n0 = n; - while (n0 > 0) { - - // retrieve shift that was temporarily stored as a negative off-diagonal element - sigma = (n0 == n) ? 0 : -work[4 * n0 - 2]; - sigmaLow = 0; - - // find start of a new split segment to process - double offDiagMin = (i0 == n0) ? 0 : work[4 * n0 - 6]; - double offDiagMax = 0; - double diagMax = work[4 * n0 - 4]; - double diagMin = diagMax; - i0 = 0; - for (int i = 4 * (n0 - 2); i >= 0; i -= 4) { - if (work[i + 2] <= 0) { - i0 = 1 + i / 4; - break; - } - if (diagMin >= 4 * offDiagMax) { - diagMin = Math.min(diagMin, work[i + 4]); - offDiagMax = Math.max(offDiagMax, work[i + 2]); - } - diagMax = Math.max(diagMax, work[i] + work[i + 2]); - offDiagMin = Math.min(offDiagMin, work[i + 2]); - } - work[4 * n0 - 2] = offDiagMin; - - // lower bound of Gershgorin disk - dMin = -Math.max(0, diagMin - 2 * Math.sqrt(diagMin * offDiagMax)); - - pingPong = 0; - int maxIter = 30 * (n0 - i0); - for (int k = 0; i0 < n0; ++k) { - if (k >= maxIter) { - throw new InvalidMatrixException(new MaxIterationsExceededException(maxIter)); - } - - // perform one step - n0 = goodStep(i0, n0); - pingPong = 1 - pingPong; - - // check for new splits after "ping" steps - // when the last elements of qd array are very small - if ((pingPong == 0) && (n0 - i0 > 3) && - (work[4 * n0 - 1] <= TOLERANCE_2 * diagMax) && - (work[4 * n0 - 2] <= TOLERANCE_2 * sigma)) { - int split = i0 - 1; - diagMax = work[4 * i0]; - offDiagMin = work[4 * i0 + 2]; - double previousEMin = work[4 * i0 + 3]; - for (int i = 4 * i0; i < 4 * n0 - 16; i += 4) { - if ((work[i + 3] <= TOLERANCE_2 * work[i]) || - (work[i + 2] <= TOLERANCE_2 * sigma)) { - // insert a split - work[i + 2] = -sigma; - split = i / 4; - diagMax = 0; - offDiagMin = work[i + 6]; - previousEMin = work[i + 7]; + if (m != j) { + if (its == maxIter) + throw new InvalidMatrixException( + new MaxIterationsExceededException(maxIter)); + its++; + double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]); + double t = Math.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 (Math.abs(p) >= Math.abs(q)) { + c = q / p; + t = Math.sqrt(c * c + 1.0); + e[i + 1] = p * t; + s = 1.0 / t; + c = c * s; } else { - diagMax = Math.max(diagMax, work[i + 4]); - offDiagMin = Math.min(offDiagMin, work[i + 2]); - previousEMin = Math.min(previousEMin, work[i + 3]); + s = p / q; + t = Math.sqrt(s * s + 1.0); + e[i + 1] = q * t; + c = 1.0 / t; + s = s * c; } - } - work[4 * n0 - 2] = offDiagMin; - work[4 * n0 - 1] = previousEMin; - i0 = split + 1; - } - } - - } - - } - - /** - * Perform two iterations with Li's tests for initial splits. - * @param n number of rows of the matrix to process - */ - private void initialSplits(final int n) { - - pingPong = 0; - for (int k = 0; k < 2; ++k) { - - // apply Li's reverse test - double d = work[4 * (n - 1) + pingPong]; - for (int i = 4 * (n - 2) + pingPong; i >= 0; i -= 4) { - if (work[i + 2] <= TOLERANCE_2 * d) { - work[i + 2] = -0.0; - d = work[i]; - } else { - d *= work[i] / (d + work[i + 2]); - } - } - - // apply dqd plus Li's forward test. - d = work[pingPong]; - for (int i = 2 + pingPong; i < 4 * n - 2; i += 4) { - final int j = i - 2 * pingPong - 1; - work[j] = d + work[i]; - if (work[i] <= TOLERANCE_2 * d) { - work[i] = -0.0; - work[j] = d; - work[j + 2] = 0.0; - d = work[i + 2]; - } else if ((MathUtils.SAFE_MIN * work[i + 2] < work[j]) && - (MathUtils.SAFE_MIN * work[j] < work[i + 2])) { - final double tmp = work[i + 2] / work[j]; - work[j + 2] = work[i] * tmp; - d *= tmp; - } else { - work[j + 2] = work[i + 2] * (work[i] / work[j]); - d *= work[i + 2] / work[j]; - } - } - work[4 * n - 3 - pingPong] = d; - - // from ping to pong - pingPong = 1 - pingPong; - - } - - } - - /** - * Perform one "good" dqd/dqds step. - *This implementation is based on Beresford N. Parlett - * and Osni A. Marques paper An - * Implementation of the dqds Algorithm (Positive Case) and on the - * corresponding LAPACK routine DLAZQ3.
- * @param start start index - * @param end end index - * @return new end (maybe deflated) - */ - private int goodStep(final int start, final int end) { - - g = 0.0; - - // step 1: accepting realEigenvalues - int deflatedEnd = end; - for (boolean deflating = true; deflating;) { - - if (start >= deflatedEnd) { - // the array has been completely deflated - return deflatedEnd; - } - - final int k = 4 * deflatedEnd + pingPong - 1; - - if ((start == deflatedEnd - 1) || - ((start != deflatedEnd - 2) && - ((work[k - 5] <= TOLERANCE_2 * (sigma + work[k - 3])) || - (work[k - 2 * pingPong - 4] <= TOLERANCE_2 * work[k - 7])))) { - - // one eigenvalue found, deflate array - work[4 * deflatedEnd - 4] = sigma + work[4 * deflatedEnd - 4 + pingPong]; - deflatedEnd -= 1; - - } else if ((start == deflatedEnd - 2) || - (work[k - 9] <= TOLERANCE_2 * sigma) || - (work[k - 2 * pingPong - 8] <= TOLERANCE_2 * work[k - 11])) { - - // two realEigenvalues found, deflate array - if (work[k - 3] > work[k - 7]) { - final double tmp = work[k - 3]; - work[k - 3] = work[k - 7]; - work[k - 7] = tmp; - } - - if (work[k - 5] > TOLERANCE_2 * work[k - 3]) { - double t = 0.5 * ((work[k - 7] - work[k - 3]) + work[k - 5]); - double s = work[k - 3] * (work[k - 5] / t); - if (s <= t) { - s = work[k - 3] * work[k - 5] / (t * (1 + Math.sqrt(1 + s / t))); - } else { - s = work[k - 3] * work[k - 5] / (t + Math.sqrt(t * (t + s))); - } - t = work[k - 7] + (s + work[k - 5]); - work[k - 3] *= work[k - 7] / t; - work[k - 7] = t; - } - work[4 * deflatedEnd - 8] = sigma + work[k - 7]; - work[4 * deflatedEnd - 4] = sigma + work[k - 3]; - deflatedEnd -= 2; - } else { - - // no more realEigenvalues found, we need to iterate - deflating = false; - - } - - } - - final int l = 4 * deflatedEnd + pingPong - 1; - - // step 2: flip array if needed - if ((dMin <= 0) || (deflatedEnd < end)) { - if (flipAllIfWarranted(deflatedEnd)) { - dMin2 = Math.min(dMin2, work[l - 1]); - work[l - 1] = - Math.min(work[l - 1], - Math.min(work[3 + pingPong], work[7 + pingPong])); - work[l - 2 * pingPong] = - Math.min(work[l - 2 * pingPong], - Math.min(work[6 + pingPong], work[6 + pingPong])); - qMax = Math.max(qMax, Math.max(work[3 + pingPong], work[7 + pingPong])); - dMin = -0.0; - } - } - - if ((dMin < 0) || - (MathUtils.SAFE_MIN * qMax < Math.min(work[l - 1], - Math.min(work[l - 9], - dMin2 + work[l - 2 * pingPong])))) { - // step 3: choose a shift - computeShiftIncrement(start, deflatedEnd, end - deflatedEnd); - - // step 4a: dqds - for (boolean loop = true; loop;) { - - // perform one dqds step with the chosen shift - dqds(start, deflatedEnd); - - // check result of the dqds step - if ((dMin >= 0) && (dMin1 > 0)) { - // the shift was good - updateSigma(tau); - return deflatedEnd; - } else if ((dMin < 0.0) && - (dMin1 > 0.0) && - (work[4 * deflatedEnd - 5 - pingPong] < TOLERANCE * (sigma + dN1)) && - (Math.abs(dN) < TOLERANCE * sigma)) { - // convergence hidden by negative DN. - work[4 * deflatedEnd - 3 - pingPong] = 0.0; - dMin = 0.0; - updateSigma(tau); - return deflatedEnd; - } else if (dMin < 0.0) { - // tau too big. Select new tau and try again. - if (tType < -22) { - // failed twice. Play it safe. - tau = 0.0; - } else if (dMin1 > 0.0) { - // late failure. Gives excellent shift. - tau = (tau + dMin) * (1.0 - 2.0 * MathUtils.EPSILON); - tType -= 11; - } else { - // early failure. Divide by 4. - tau *= 0.25; - tType -= 12; - } - } else if (Double.isNaN(dMin)) { - tau = 0.0; - } else { - // possible underflow. Play it safe. - loop = false; - } - } - - } - - // perform a dqd step (i.e. no shift) - dqd(start, deflatedEnd); - - return deflatedEnd; - - } - - /** - * Flip all elements of qd array if warranted. - * @param n number of rows in the block - * @return true if qd array was flipped - */ - private boolean flipAllIfWarranted(final int n) { - if (1.5 * work[pingPong] >= work[4 * (n - 1) + pingPong]) { - return false; - } - - int j = 4 * (n - 1); - for (int i = 0; i < j; i += 4) { - final double tmp1 = work[i]; - work[i] = work[j]; - work[j] = tmp1; - final double tmp2 = work[i+1]; - work[i+1] = work[j+1]; - work[j+1] = tmp2; - final double tmp3 = work[i+2]; - work[i+2] = work[j-2]; - work[j-2] = tmp3; - final double tmp4 = work[i+3]; - work[i+3] = work[j-1]; - work[j-1] = tmp4; - j -= 4; - } - - return true; - - } - - /** - * Flip every other elements of qd array if warranted. - * @param n number of rows in the block - * @return true if qd array was flipped - */ - private boolean flipEveryOtherIfWarranted(final int n) { - if (1.5 * work[pingPong] >= work[4 * (n - 1) + pingPong]) { - return false; - } - - // flip array - int j = 4 * (n - 1); - for (int i = 0; i < j; i += 4) { - for (int k = 0; k < 4; k += 2) { - final double tmp = work[i + k]; - work[i + k] = work[j - k]; - work[j - k] = tmp; - } - j -= 4; - } - - return true; - - } - - /** - * Compute an interval containing all realEigenvalues of a block. - * @param index index of the first row of the block - * @param n number of rows of the block - * @return an interval containing the realEigenvalues - */ - private double[] eigenvaluesRange(final int index, final int n) { - - // find the bounds of the spectra of the local block - final int lowerStart = 4 * main.length; - final int upperStart = 5 * main.length; - double lower = Double.POSITIVE_INFINITY; - double upper = Double.NEGATIVE_INFINITY; - for (int i = 0; i < n; ++i) { - lower = Math.min(lower, work[lowerStart + index +i]); - upper = Math.max(upper, work[upperStart + index +i]); - } - - // set thresholds - final double tNorm = Math.max(Math.abs(lower), Math.abs(upper)); - final double relativeTolerance = Math.sqrt(MathUtils.EPSILON); - final double absoluteTolerance = 4 * minPivot; - final int maxIter = - 2 + (int) ((Math.log(tNorm + minPivot) - Math.log(minPivot)) / Math.log(2.0)); - final double margin = 2 * (tNorm * MathUtils.EPSILON * n + 2 * minPivot); - - // search lower eigenvalue - double left = lower - margin; - double right = upper + margin; - for (int i = 0; i < maxIter; ++i) { - - final double range = right - left; - if ((range < absoluteTolerance) || - (range < relativeTolerance * Math.max(Math.abs(left), Math.abs(right)))) { - // search has converged - break; - } - - final double middle = 0.5 * (left + right); - if (countEigenValues(middle, index, n) >= 1) { - right = middle; - } else { - left = middle; - } - - } - lower = Math.max(lower, left - 100 * MathUtils.EPSILON * Math.abs(left)); - - // search upper eigenvalue - left = lower - margin; - right = upper + margin; - for (int i = 0; i < maxIter; ++i) { - - final double range = right - left; - if ((range < absoluteTolerance) || - (range < relativeTolerance * Math.max(Math.abs(left), Math.abs(right)))) { - // search has converged - break; - } - - final double middle = 0.5 * (left + right); - if (countEigenValues(middle, index, n) >= n) { - right = middle; - } else { - left = middle; - } - - } - upper = Math.min(upper, right + 100 * MathUtils.EPSILON * Math.abs(right)); - - return new double[] { lower, upper }; - - } - - /** - * Count the number of realEigenvalues below a point. - * @param t value below which we must count the number of realEigenvalues - * @param index index of the first row of the block - * @param n number of rows of the block - * @return number of realEigenvalues smaller than t - */ - private int countEigenValues(final double t, final int index, final int n) { - double ratio = main[index] - t; - int count = (ratio > 0) ? 0 : 1; - for (int i = 1; i < n; ++i) { - ratio = main[index + i] - squaredSecondary[index + i - 1] / ratio - t; - if (ratio <= 0) { - ++count; - } - } - return count; - } - - /** - * Decompose the shifted tridiagonal matrix T-λI as LDLT. - *A shifted symmetric tridiagonal matrix T can be decomposed as - * LDLT where L is a lower bidiagonal matrix with unit diagonal - * and D is a diagonal matrix. This method is an implementation of - * algorithm 4.4.7 from Dhillon's thesis.
- * @param lambda shift to add to the matrix before decomposing it - * to ensure it is positive definite - * @param index index of the first row of the block - * @param n number of rows of the block - */ - private void ldlTDecomposition(final double lambda, final int index, final int n) { - double di = main[index] - lambda; - work[0] = Math.abs(di); - for (int i = 1; i < n; ++i) { - final int fourI = 4 * i; - final double eiM1 = secondary[index + i - 1]; - final double ratio = eiM1 / di; - work[fourI - 2] = ratio * ratio * Math.abs(di); - di = (main[index + i] - lambda) - eiM1 * ratio; - work[fourI] = Math.abs(di); - } - } - - /** - * Perform a dqds step, using current shift increment. - *This implementation is a translation of the LAPACK routine DLASQ5.
- * @param start start index - * @param end end index - */ - private void dqds(final int start, final int end) { - - eMin = work[4 * start + pingPong + 4]; - double d = work[4 * start + pingPong] - tau; - dMin = d; - dMin1 = -work[4 * start + pingPong]; - - if (pingPong == 0) { - for (int j4 = 4 * start + 3; j4 <= 4 * (end - 3); j4 += 4) { - work[j4 - 2] = d + work[j4 - 1]; - final double tmp = work[j4 + 1] / work[j4 - 2]; - d = d * tmp - tau; - dMin = Math.min(dMin, d); - work[j4] = work[j4 - 1] * tmp; - eMin = Math.min(work[j4], eMin); - } - } else { - for (int j4 = 4 * start + 3; j4 <= 4 * (end - 3); j4 += 4) { - work[j4 - 3] = d + work[j4]; - final double tmp = work[j4 + 2] / work[j4 - 3]; - d = d * tmp - tau; - dMin = Math.min(dMin, d); - work[j4 - 1] = work[j4] * tmp; - eMin = Math.min(work[j4 - 1], eMin); - } - } - - // unroll last two steps. - dN2 = d; - dMin2 = dMin; - int j4 = 4 * (end - 2) - pingPong - 1; - int j4p2 = j4 + 2 * pingPong - 1; - work[j4 - 2] = dN2 + work[j4p2]; - work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]); - dN1 = work[j4p2 + 2] * (dN2 / work[j4 - 2]) - tau; - dMin = Math.min(dMin, dN1); - - dMin1 = dMin; - j4 = j4 + 4; - j4p2 = j4 + 2 * pingPong - 1; - work[j4 - 2] = dN1 + work[j4p2]; - work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]); - dN = work[j4p2 + 2] * (dN1 / work[j4 - 2]) - tau; - dMin = Math.min(dMin, dN); - - work[j4 + 2] = dN; - work[4 * end - pingPong - 1] = eMin; - - } - - - /** - * Perform a dqd step. - *This implementation is a translation of the LAPACK routine DLASQ6.
- * @param start start index - * @param end end index - */ - private void dqd(final int start, final int end) { - - eMin = work[4 * start + pingPong + 4]; - double d = work[4 * start + pingPong]; - dMin = d; - - if (pingPong == 0) { - for (int j4 = 4 * start + 3; j4 < 4 * (end - 3); j4 += 4) { - work[j4 - 2] = d + work[j4 - 1]; - if (work[j4 - 2] == 0.0) { - work[j4] = 0.0; - d = work[j4 + 1]; - dMin = d; - eMin = 0.0; - } else if ((MathUtils.SAFE_MIN * work[j4 + 1] < work[j4 - 2]) && - (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4 + 1])) { - final double tmp = work[j4 + 1] / work[j4 - 2]; - work[j4] = work[j4 - 1] * tmp; - d *= tmp; - } else { - work[j4] = work[j4 + 1] * (work[j4 - 1] / work[j4 - 2]); - d *= work[j4 + 1] / work[j4 - 2]; - } - dMin = Math.min(dMin, d); - eMin = Math.min(eMin, work[j4]); - } - } else { - for (int j4 = 4 * start + 3; j4 < 4 * (end - 3); j4 += 4) { - work[j4 - 3] = d + work[j4]; - if (work[j4 - 3] == 0.0) { - work[j4 - 1] = 0.0; - d = work[j4 + 2]; - dMin = d; - eMin = 0.0; - } else if ((MathUtils.SAFE_MIN * work[j4 + 2] < work[j4 - 3]) && - (MathUtils.SAFE_MIN * work[j4 - 3] < work[j4 + 2])) { - final double tmp = work[j4 + 2] / work[j4 - 3]; - work[j4 - 1] = work[j4] * tmp; - d *= tmp; - } else { - work[j4 - 1] = work[j4 + 2] * (work[j4] / work[j4 - 3]); - d *= work[j4 + 2] / work[j4 - 3]; - } - dMin = Math.min(dMin, d); - eMin = Math.min(eMin, work[j4 - 1]); - } - } - - // Unroll last two steps - dN2 = d; - dMin2 = dMin; - int j4 = 4 * (end - 2) - pingPong - 1; - int j4p2 = j4 + 2 * pingPong - 1; - work[j4 - 2] = dN2 + work[j4p2]; - if (work[j4 - 2] == 0.0) { - work[j4] = 0.0; - dN1 = work[j4p2 + 2]; - dMin = dN1; - eMin = 0.0; - } else if ((MathUtils.SAFE_MIN * work[j4p2 + 2] < work[j4 - 2]) && - (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4p2 + 2])) { - final double tmp = work[j4p2 + 2] / work[j4 - 2]; - work[j4] = work[j4p2] * tmp; - dN1 = dN2 * tmp; - } else { - work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]); - dN1 = work[j4p2 + 2] * (dN2 / work[j4 - 2]); - } - dMin = Math.min(dMin, dN1); - - dMin1 = dMin; - j4 = j4 + 4; - j4p2 = j4 + 2 * pingPong - 1; - work[j4 - 2] = dN1 + work[j4p2]; - if (work[j4 - 2] == 0.0) { - work[j4] = 0.0; - dN = work[j4p2 + 2]; - dMin = dN; - eMin = 0.0; - } else if ((MathUtils.SAFE_MIN * work[j4p2 + 2] < work[j4 - 2]) && - (MathUtils.SAFE_MIN * work[j4 - 2] < work[j4p2 + 2])) { - final double tmp = work[j4p2 + 2] / work[j4 - 2]; - work[j4] = work[j4p2] * tmp; - dN = dN1 * tmp; - } else { - work[j4] = work[j4p2 + 2] * (work[j4p2] / work[j4 - 2]); - dN = work[j4p2 + 2] * (dN1 / work[j4 - 2]); - } - dMin = Math.min(dMin, dN); - - work[j4 + 2] = dN; - work[4 * end - pingPong - 1] = eMin; - - } - - /** - * Compute the shift increment as an estimate of the smallest eigenvalue. - *This implementation is a translation of the LAPACK routine DLAZQ4.
- * @param start start index - * @param end end index - * @param deflated number of realEigenvalues just deflated - */ - private void computeShiftIncrement(final int start, final int end, final int deflated) { - - final double cnst1 = 0.563; - final double cnst2 = 1.010; - final double cnst3 = 1.05; - - // a negative dMin forces the shift to take that absolute value - // tType records the type of shift. - if (dMin <= 0.0) { - tau = -dMin; - tType = -1; - return; - } - - int nn = 4 * end + pingPong - 1; - switch (deflated) { - - case 0 : // no realEigenvalues deflated. - if (dMin == dN || dMin == dN1) { - - double b1 = Math.sqrt(work[nn - 3]) * Math.sqrt(work[nn - 5]); - double b2 = Math.sqrt(work[nn - 7]) * Math.sqrt(work[nn - 9]); - double a2 = work[nn - 7] + work[nn - 5]; - - if (dMin == dN && dMin1 == dN1) { - // cases 2 and 3. - final double gap2 = dMin2 - a2 - dMin2 * 0.25; - final double gap1 = a2 - dN - ((gap2 > 0.0 && gap2 > b2) ? (b2 / gap2) * b2 : (b1 + b2)); - if (gap1 > 0.0 && gap1 > b1) { - tau = Math.max(dN - (b1 / gap1) * b1, 0.5 * dMin); - tType = -2; - } else { - double s = 0.0; - if (dN > b1) { - s = dN - b1; - } - if (a2 > (b1 + b2)) { - s = Math.min(s, a2 - (b1 + b2)); - } - tau = Math.max(s, 0.333 * dMin); - tType = -3; - } - } else { - // case 4. - tType = -4; - double s = 0.25 * dMin; - double gam; - int np; - if (dMin == dN) { - gam = dN; - a2 = 0.0; - if (work[nn - 5] > work[nn - 7]) { - return; - } - b2 = work[nn - 5] / work[nn - 7]; - np = nn - 9; - } else { - np = nn - 2 * pingPong; - b2 = work[np - 2]; - gam = dN1; - if (work[np - 4] > work[np - 2]) { - return; - } - a2 = work[np - 4] / work[np - 2]; - if (work[nn - 9] > work[nn - 11]) { - return; - } - b2 = work[nn - 9] / work[nn - 11]; - np = nn - 13; - } - - // approximate contribution to norm squared from i < nn-1. - a2 = a2 + b2; - for (int i4 = np; i4 >= 4 * start + 2 + pingPong; i4 -= 4) { - if(b2 == 0.0) { + if (e[i + 1] == 0.0) { + realEigenvalues[i + 1] -= u; + e[m] = 0.0; break; } - b1 = b2; - if (work[i4] > work[i4 - 2]) { - return; - } - b2 = b2 * (work[i4] / work[i4 - 2]); - a2 = a2 + b2; - if (100 * Math.max(b2, b1) < a2 || cnst1 < a2) { - 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; } } - a2 = cnst3 * a2; - - // rayleigh quotient residual bound. - if (a2 < cnst1) { - s = gam * (1 - Math.sqrt(a2)) / (1 + a2); - } - tau = s; - + if (e[i + 1] == 0.0 && i >= j) + continue; + realEigenvalues[j] -= u; + e[j] = q; + e[m] = 0.0; } - } else if (dMin == dN2) { + } while (m != j); + } - // case 5. - tType = -5; - double s = 0.25 * dMin; - - // compute contribution to norm squared from i > nn-2. - final int np = nn - 2 * pingPong; - double b1 = work[np - 2]; - double b2 = work[np - 6]; - final double gam = dN2; - if (work[np - 8] > b2 || work[np - 4] > b1) { - return; + //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]; } - double a2 = (work[np - 8] / b2) * (1 + work[np - 4] / b1); - - // approximate contribution to norm squared from i < nn-2. - if (end - start > 3) { - b2 = work[nn - 13] / work[nn - 15]; - a2 = a2 + b2; - for (int i4 = nn - 17; i4 >= 4 * start + 2 + pingPong; i4 -= 4) { - if (b2 == 0.0) { - break; - } - b1 = b2; - if (work[i4] > work[i4 - 2]) { - return; - } - b2 = b2 * (work[i4] / work[i4 - 2]); - a2 = a2 + b2; - if (100 * Math.max(b2, b1) < a2 || cnst1 < a2) { - break; - } - } - a2 = cnst3 * a2; - } - - if (a2 < cnst1) { - tau = gam * (1 - Math.sqrt(a2)) / (1 + a2); - } else { - tau = s; - } - - } else { - - // case 6, no information to guide us. - if (tType == -6) { - g += 0.333 * (1 - g); - } else if (tType == -18) { - g = 0.25 * 0.333; - } else { - g = 0.25; - } - tau = g * dMin; - tType = -6; - } - break; - - case 1 : // one eigenvalue just deflated. use dMin1, dN1 for dMin and dN. - if (dMin1 == dN1 && dMin2 == dN2) { - - // cases 7 and 8. - tType = -7; - double s = 0.333 * dMin1; - if (work[nn - 5] > work[nn - 7]) { - return; + 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; } - double b1 = work[nn - 5] / work[nn - 7]; - double b2 = b1; - if (b2 != 0.0) { - for (int i4 = 4 * end - 10 + pingPong; i4 >= 4 * start + 2 + pingPong; i4 -= 4) { - final double oldB1 = b1; - if (work[i4] > work[i4 - 2]) { - return; - } - b1 = b1 * (work[i4] / work[i4 - 2]); - b2 = b2 + b1; - if (100 * Math.max(b1, oldB1) < b2) { - break; - } - } - } - b2 = Math.sqrt(cnst3 * b2); - final double a2 = dMin1 / (1 + b2 * b2); - final double gap2 = 0.5 * dMin2 - a2; - if (gap2 > 0.0 && gap2 > b2 * a2) { - tau = Math.max(s, a2 * (1 - cnst2 * a2 * (b2 / gap2) * b2)); - } else { - tau = Math.max(s, a2 * (1 - cnst2 * b2)); - tType = -8; - } - } else { - - // case 9. - tau = 0.25 * dMin1; - if (dMin1 == dN1) { - tau = 0.5 * dMin1; - } - tType = -9; - } - break; - - case 2 : // two realEigenvalues deflated. use dMin2, dN2 for dMin and dN. - - // cases 10 and 11. - if (dMin2 == dN2 && 2 * work[nn - 5] < work[nn - 7]) { - tType = -10; - final double s = 0.333 * dMin2; - if (work[nn - 5] > work[nn - 7]) { - return; - } - double b1 = work[nn - 5] / work[nn - 7]; - double b2 = b1; - if (b2 != 0.0){ - for (int i4 = 4 * end - 9 + pingPong; i4 >= 4 * start + 2 + pingPong; i4 -= 4) { - if (work[i4] > work[i4 - 2]) { - return; - } - b1 *= work[i4] / work[i4 - 2]; - b2 += b1; - if (100 * b1 < b2) { - break; - } - } - } - b2 = Math.sqrt(cnst3 * b2); - final double a2 = dMin2 / (1 + b2 * b2); - final double gap2 = work[nn - 7] + work[nn - 9] - - Math.sqrt(work[nn - 11]) * Math.sqrt(work[nn - 9]) - a2; - if (gap2 > 0.0 && gap2 > b2 * a2) { - tau = Math.max(s, a2 * (1 - cnst2 * a2 * (b2 / gap2) * b2)); - } else { - tau = Math.max(s, a2 * (1 - cnst2 * b2)); - } - } else { - tau = 0.25 * dMin2; - tType = -11; - } - break; - - default : // case 12, more than two realEigenvalues deflated. no information. - tau = 0.0; - tType = -12; - } - - } - - /** - * Update sigma. - * @param shift shift to apply to sigma - */ - private void updateSigma(final double shift) { - // BEWARE: do NOT attempt to simplify the following statements - // the expressions below take care to accumulate the part of sigma - // that does not fit within a double variable into sigmaLow - if (shift < sigma) { - sigmaLow += shift; - final double t = sigma + sigmaLow; - sigmaLow -= t - sigma; - sigma = t; - } else { - final double t = sigma + shift; - sigmaLow += sigma - (t - shift); - sigma = t; - } - } - - /** - * Find eigenvectors. - */ - private void findEigenVectors() { - - final int m = main.length; - eigenvectors = new ArrayRealVector[m]; - - // perform an initial non-shifted LDLt decomposition - final double[] d = new double[m]; - final double[] l = new double[m - 1]; - // avoid zero divide on indefinite matrix - final double mu = realEigenvalues[m-1] <= 0 && realEigenvalues[0] > 0 ? 0.5-realEigenvalues[m-1] : 0; - double di = main[0]+mu; - d[0] = di; - for (int i = 1; i < m; ++i) { - final double eiM1 = secondary[i - 1]; - final double ratio = eiM1 / di; - di = main[i] - eiM1 * ratio + mu; - l[i - 1] = ratio; - d[i] = di; - } - - // compute eigenvectors - for (int i = 0; i < m; ++i) { - eigenvectors[i] = findEigenvector(realEigenvalues[i]+mu, d, l); - } - - } - - /** - * Find an eigenvector corresponding to an eigenvalue, using bidiagonals. - *This method corresponds to algorithm X from Dhillon's thesis.
- * - * @param eigenvalue eigenvalue for which eigenvector is desired - * @param d diagonal elements of the initial non-shifted D matrix - * @param l off-diagonal elements of the initial non-shifted L matrix - * @return an eigenvector - */ - private ArrayRealVector findEigenvector(final double eigenvalue, - final double[] d, final double[] l) { - - // compute the LDLt and UDUt decompositions of the - // perfectly shifted tridiagonal matrix - final int m = main.length; - stationaryQuotientDifferenceWithShift(d, l, eigenvalue); - progressiveQuotientDifferenceWithShift(d, l, eigenvalue); - - // select the twist index leading to - // the least diagonal element in the twisted factorization - int r = m - 1; - double minG = Math.abs(work[6 * r] + work[6 * r + 3] + eigenvalue); - int sixI = 0; - for (int i = 0; i < m - 1; ++i) { - final double absG = Math.abs(work[sixI] + d[i] * work[sixI + 9] / work[sixI + 10]); - if (absG < minG) { - r = i; - minG = absG; - } - sixI += 6; - } - - // solve the singular system by ignoring the equation - // at twist index and propagating upwards and downwards - double[] eigenvector = new double[m]; - double n2 = 1; - eigenvector[r] = 1; - double z = 1; - for (int i = r - 1; i >= 0; --i) { - z *= -work[6 * i + 2]; - eigenvector[i] = z; - n2 += z * z; - } - z = 1; - for (int i = r + 1; i < m; ++i) { - z *= -work[6 * i - 1]; - eigenvector[i] = z; - n2 += z * z; - } - - // normalize vector - final double inv = 1.0 / Math.sqrt(n2); - for (int i = 0; i < m; ++i) { - eigenvector[i] *= inv; - } - - return (transformer == null) ? - new ArrayRealVector(eigenvector, false) : - new ArrayRealVector(transformer.getQ().operate(eigenvector), false); - - } - - /** - * Decompose matrix LDLT - λ I as - * L+D+L+T. - *This method corresponds to algorithm 4.4.3 (dstqds) from Dhillon's thesis.
- * @param d diagonal elements of D, - * @param l off-diagonal elements of L - * @param lambda shift to apply - */ - private void stationaryQuotientDifferenceWithShift(final double[] d, final double[] l, - final double lambda) { - final int nM1 = d.length - 1; - double si = -lambda; - int sixI = 0; - for (int i = 0; i < nM1; ++i) { - final double di = d[i]; - final double li = l[i]; - final double ldi = li * di; - final double diP1 = di + si; - final double liP1 = ldi / diP1; - work[sixI] = si; - work[sixI + 1] = diP1; - work[sixI + 2] = liP1; - si = li * liP1 * si - lambda; - sixI += 6; - } - if (Double.isNaN(si)) { - // one of the pivot was null, use a slower but safer version of dstqds - si = -lambda; - sixI = 0; - for (int i = 0; i < nM1; ++i) { - final double di = d[i]; - final double li = l[i]; - final double ldi = li * di; - double diP1 = di + si; - if (Math.abs(diP1) < minPivot) { - diP1 = -minPivot; - } - final double liP1 = ldi / diP1; - work[sixI] = si; - work[sixI + 1] = diP1; - work[sixI + 2] = liP1; - si = li * ((liP1 == 0) ? li * di : liP1 * si) - lambda; - sixI += 6; } } - work[6 * nM1 + 1] = d[nM1] + si; - work[6 * nM1] = si; - } - /** - * Decompose matrix LDLT - λ I as - * U-D-U-T. - *This method corresponds to algorithm 4.4.5 (dqds) from Dhillon's thesis.
- * @param d diagonal elements of D - * @param l off-diagonal elements of L - * @param lambda shift to apply - */ - private void progressiveQuotientDifferenceWithShift(final double[] d, final double[] l, - final double lambda) { - final int nM1 = d.length - 1; - double pi = d[nM1] - lambda; - int sixI = 6 * (nM1 - 1); - for (int i = nM1 - 1; i >= 0; --i) { - final double di = d[i]; - final double li = l[i]; - final double diP1 = di * li * li + pi; - final double t = di / diP1; - work[sixI + 9] = pi; - work[sixI + 10] = diP1; - work[sixI + 5] = li * t; - pi = pi * t - lambda; - sixI -= 6; - } - if (Double.isNaN(pi)) { - // one of the pivot was null, use a slower but safer version of dqds - pi = d[nM1] - lambda; - sixI = 6 * (nM1 - 1); - for (int i = nM1 - 1; i >= 0; --i) { - final double di = d[i]; - final double li = l[i]; - double diP1 = di * li * li + pi; - if (Math.abs(diP1) < minPivot) { - diP1 = -minPivot; - } - final double t = di / diP1; - work[sixI + 9] = pi; - work[sixI + 10] = diP1; - work[sixI + 5] = li * t; - pi = ((t == 0) ? di : pi * t) - lambda; - sixI -= 6; + // Determine the largest eigen value in absolute term. + double maxAbsoluteValue=0.0; + for (int i = 0; i < n; i++) { + if (Math.abs(realEigenvalues[i])>maxAbsoluteValue) { + maxAbsoluteValue=Math.abs(realEigenvalues[i]); } } - work[3] = pi; - work[4] = pi; + // Make null any eigen value too small to be significant + if (maxAbsoluteValue!=0.0) { + for (int i=0; i < n; i++) { + if (Math.abs(realEigenvalues[i])- * Note that since this class computes only the compact or truncated SVD and not - * the full SVD, the singular values computed are always positive. + * The Singular Value Decomposition of matrix A is a set of three matrices: U, + * Σ and V such that A = U × Σ × VT. Let A be + * a m × n matrix, then U is a m × p orthogonal matrix, Σ is a + * p × p diagonal matrix with positive or null elements, V is a p × + * n orthogonal matrix (hence VT is also orthogonal) where + * p=min(m,n). *
- * * @version $Revision$ $Date$ * @since 2.0 */ -public class SingularValueDecompositionImpl implements SingularValueDecomposition { +public class SingularValueDecompositionImpl implements + SingularValueDecomposition { /** Number of rows of the initial matrix. */ private int m; @@ -51,21 +41,6 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio /** Number of columns of the initial matrix. */ private int n; - /** Transformer to bidiagonal. */ - private BiDiagonalTransformer transformer; - - /** Main diagonal of the bidiagonal matrix. */ - private double[] mainBidiagonal; - - /** Secondary diagonal of the bidiagonal matrix. */ - private double[] secondaryBidiagonal; - - /** Main diagonal of the tridiagonal matrix. */ - private double[] mainTridiagonal; - - /** Secondary diagonal of the tridiagonal matrix. */ - private double[] secondaryTridiagonal; - /** Eigen decomposition of the tridiagonal matrix. */ private EigenDecomposition eigenDecomposition; @@ -89,120 +64,101 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio /** * Calculates the compact Singular Value Decomposition of the given matrix. - * @param matrix The matrix to decompose. - * @exception InvalidMatrixException (wrapping a {@link - * org.apache.commons.math.ConvergenceException} if algorithm fails to converge + * @param matrix + * The matrix to decompose. + * @exception InvalidMatrixException + * (wrapping a + * {@link org.apache.commons.math.ConvergenceException} if + * algorithm fails to converge */ public SingularValueDecompositionImpl(final RealMatrix matrix) - throws InvalidMatrixException { - this(matrix, Math.min(matrix.getRowDimension(), matrix.getColumnDimension())); - } - - /** - * Calculates the Singular Value Decomposition of the given matrix. - * @param matrix The matrix to decompose. - * @param max maximal number of singular values to compute - * @exception InvalidMatrixException (wrapping a {@link - * org.apache.commons.math.ConvergenceException} if algorithm fails to converge - */ - public SingularValueDecompositionImpl(final RealMatrix matrix, final int max) - throws InvalidMatrixException { + throws InvalidMatrixException { m = matrix.getRowDimension(); n = matrix.getColumnDimension(); - cachedU = null; - cachedS = null; - cachedV = null; + cachedU = null; + cachedS = null; + cachedV = null; cachedVt = null; - // transform the matrix to bidiagonal - transformer = new BiDiagonalTransformer(matrix); - mainBidiagonal = transformer.getMainDiagonalRef(); - secondaryBidiagonal = transformer.getSecondaryDiagonalRef(); - - // compute Bt.B (if upper diagonal) or B.Bt (if lower diagonal) - mainTridiagonal = new double[mainBidiagonal.length]; - secondaryTridiagonal = new double[mainBidiagonal.length - 1]; - double a = mainBidiagonal[0]; - mainTridiagonal[0] = a * a; - for (int i = 1; i < mainBidiagonal.length; ++i) { - final double b = secondaryBidiagonal[i - 1]; - secondaryTridiagonal[i - 1] = a * b; - a = mainBidiagonal[i]; - mainTridiagonal[i] = a * a + b * b; + double[][] localcopy = matrix.getData(); + double[][] matATA = new double[n][n]; + // + // create A^T*A + // + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + matATA[i][j] = 0.0; + for (int k = 0; k < m; k++) { + matATA[i][j] += localcopy[k][i] * localcopy[k][j]; + } + } } - // compute singular values - eigenDecomposition = - new EigenDecompositionImpl(mainTridiagonal, secondaryTridiagonal, - MathUtils.SAFE_MIN); - final double[] eigenValues = eigenDecomposition.getRealEigenvalues(); - int p = Math.min(max, eigenValues.length); - while ((p > 0) && (eigenValues[p - 1] <= 0)) { - --p; - } - singularValues = new double[p]; - for (int i = 0; i < p; ++i) { - singularValues[i] = Math.sqrt(eigenValues[i]); + double[][] matAAT = new double[m][m]; + // + // create A*A^T + // + for (int i = 0; i < m; i++) { + for (int j = 0; j < m; j++) { + matAAT[i][j] = 0.0; + for (int k = 0; k < n; k++) { + matAAT[i][j] += localcopy[i][k] * localcopy[j][k]; + } + } } + int p; + if (m>=n) { + p=n; + // compute eigen decomposition of A^T*A + eigenDecomposition = new EigenDecompositionImpl( + new Array2DRowRealMatrix(matATA),1.0); + singularValues = eigenDecomposition.getRealEigenvalues(); + cachedV = eigenDecomposition.getV(); + // compute eigen decomposition of A*A^T + eigenDecomposition = new EigenDecompositionImpl( + new Array2DRowRealMatrix(matAAT),1.0); + cachedU = eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1); + } else { + p=m; + // compute eigen decomposition of A*A^T + eigenDecomposition = new EigenDecompositionImpl( + new Array2DRowRealMatrix(matAAT),1.0); + singularValues = eigenDecomposition.getRealEigenvalues(); + cachedU = eigenDecomposition.getV(); + + // compute eigen decomposition of A^T*A + eigenDecomposition = new EigenDecompositionImpl( + new Array2DRowRealMatrix(matATA),1.0); + cachedV = eigenDecomposition.getV().getSubMatrix(0,n-1,0,p-1); + } + for (int i = 0; i < p; i++) { + singularValues[i] = Math.sqrt(Math.abs(singularValues[i])); + } + // Up to this point, U and V are computed independently of each other. + // There still an sign indetermination of each column of, say, U. + // The sign is set such that A.V_i=sigma_i.U_i (i<=p) + // The right sign corresponds to a positive dot product of A.V_i and U_i + for (int i = 0; i < p; i++) { + RealVector tmp = cachedU.getColumnVector(i); + double product=matrix.operate(cachedV.getColumnVector(i)).dotProduct(tmp); + if (product<0) { + cachedU.setColumnVector(i, tmp.mapMultiply(-1.0)); + } + } } /** {@inheritDoc} */ - public RealMatrix getU() - throws InvalidMatrixException { - - if (cachedU == null) { - - final int p = singularValues.length; - if (m >= n) { - // the tridiagonal matrix is Bt.B, where B is upper bidiagonal - final RealMatrix e = - eigenDecomposition.getV().getSubMatrix(0, n - 1, 0, p - 1); - final double[][] eData = e.getData(); - final double[][] wData = new double[m][p]; - double[] ei1 = eData[0]; - for (int i = 0; i < p; ++i) { - // compute W = B.E.S^(-1) where E is the eigenvectors matrix - final double mi = mainBidiagonal[i]; - final double[] ei0 = ei1; - final double[] wi = wData[i]; - if (i < n - 1) { - ei1 = eData[i + 1]; - final double si = secondaryBidiagonal[i]; - for (int j = 0; j < p; ++j) { - wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j]; - } - } else { - for (int j = 0; j < p; ++j) { - wi[j] = mi * ei0[j] / singularValues[j]; - } - } - } - - for (int i = p; i < m; ++i) { - wData[i] = new double[p]; - } - cachedU = - transformer.getU().multiply(MatrixUtils.createRealMatrix(wData)); - } else { - // the tridiagonal matrix is B.Bt, where B is lower bidiagonal - final RealMatrix e = - eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1); - cachedU = transformer.getU().multiply(e); - } - - } - + public RealMatrix getU() throws InvalidMatrixException { // return the cached matrix return cachedU; } /** {@inheritDoc} */ - public RealMatrix getUT() - throws InvalidMatrixException { + public RealMatrix getUT() throws InvalidMatrixException { if (cachedUt == null) { cachedUt = getU().transpose(); @@ -214,8 +170,7 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio } /** {@inheritDoc} */ - public RealMatrix getS() - throws InvalidMatrixException { + public RealMatrix getS() throws InvalidMatrixException { if (cachedS == null) { @@ -227,64 +182,19 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio } /** {@inheritDoc} */ - public double[] getSingularValues() - throws InvalidMatrixException { + public double[] getSingularValues() throws InvalidMatrixException { return singularValues.clone(); } /** {@inheritDoc} */ - public RealMatrix getV() - throws InvalidMatrixException { - - if (cachedV == null) { - - final int p = singularValues.length; - if (m >= n) { - // the tridiagonal matrix is Bt.B, where B is upper bidiagonal - final RealMatrix e = - eigenDecomposition.getV().getSubMatrix(0, n - 1, 0, p - 1); - cachedV = transformer.getV().multiply(e); - } else { - // the tridiagonal matrix is B.Bt, where B is lower bidiagonal - // compute W = Bt.E.S^(-1) where E is the eigenvectors matrix - final RealMatrix e = - eigenDecomposition.getV().getSubMatrix(0, m - 1, 0, p - 1); - final double[][] eData = e.getData(); - final double[][] wData = new double[n][p]; - double[] ei1 = eData[0]; - for (int i = 0; i < p; ++i) { - final double mi = mainBidiagonal[i]; - final double[] ei0 = ei1; - final double[] wi = wData[i]; - if (i < m - 1) { - ei1 = eData[i + 1]; - final double si = secondaryBidiagonal[i]; - for (int j = 0; j < p; ++j) { - wi[j] = (mi * ei0[j] + si * ei1[j]) / singularValues[j]; - } - } else { - for (int j = 0; j < p; ++j) { - wi[j] = mi * ei0[j] / singularValues[j]; - } - } - } - for (int i = p; i < n; ++i) { - wData[i] = new double[p]; - } - cachedV = - transformer.getV().multiply(MatrixUtils.createRealMatrix(wData)); - } - - } - + public RealMatrix getV() throws InvalidMatrixException { // return the cached matrix return cachedV; } /** {@inheritDoc} */ - public RealMatrix getVT() - throws InvalidMatrixException { + public RealMatrix getVT() throws InvalidMatrixException { if (cachedVt == null) { cachedVt = getV().transpose(); @@ -307,15 +217,16 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio if (dimension == 0) { throw MathRuntimeException.createIllegalArgumentException( - "cutoff singular value is {0}, should be at most {1}", - minSingularValue, singularValues[0]); + "cutoff singular value is {0}, should be at most {1}", + minSingularValue, singularValues[0]); } final double[][] data = new double[dimension][p]; getVT().walkInOptimizedOrder(new DefaultRealMatrixPreservingVisitor() { /** {@inheritDoc} */ @Override - public void visit(final int row, final int column, final double value) { + public void visit(final int row, final int column, + final double value) { data[row][column] = value / singularValues[row]; } }, 0, dimension - 1, 0, p - 1); @@ -326,27 +237,24 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio } /** {@inheritDoc} */ - public double getNorm() - throws InvalidMatrixException { + public double getNorm() throws InvalidMatrixException { return singularValues[0]; } /** {@inheritDoc} */ - public double getConditionNumber() - throws InvalidMatrixException { + public double getConditionNumber() throws InvalidMatrixException { return singularValues[0] / singularValues[singularValues.length - 1]; } /** {@inheritDoc} */ - public int getRank() - throws IllegalStateException { + public int getRank() throws IllegalStateException { final double threshold = Math.max(m, n) * Math.ulp(singularValues[0]); for (int i = singularValues.length - 1; i >= 0; --i) { - if (singularValues[i] > threshold) { - return i + 1; - } + if (singularValues[i] > threshold) { + return i + 1; + } } return 0; @@ -354,8 +262,8 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio /** {@inheritDoc} */ public DecompositionSolver getSolver() { - return new Solver(singularValues, getUT(), getV(), - getRank() == Math.max(m, n)); + return new Solver(singularValues, getUT(), getV(), getRank() == Math + .max(m, n)); } /** Specialized solver. */ @@ -369,58 +277,81 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio /** * Build a solver from decomposed matrix. - * @param singularValues singularValues - * @param uT UT matrix of the decomposition - * @param v V matrix of the decomposition - * @param nonSingular singularity indicator + * @param singularValues + * singularValues + * @param uT + * UT matrix of the decomposition + * @param v + * V matrix of the decomposition + * @param nonSingular + * singularity indicator */ - private Solver(final double[] singularValues, final RealMatrix uT, final RealMatrix v, - final boolean nonSingular) { - double[][] suT = uT.getData(); + private Solver(final double[] singularValues, final RealMatrix uT, + final RealMatrix v, final boolean nonSingular) { + double[][] suT = uT.getData(); for (int i = 0; i < singularValues.length; ++i) { - final double a = 1.0 / singularValues[i]; + final double a; + if (singularValues[i]>0) { + a=1.0 / singularValues[i]; + } else { + a=0.0; + } final double[] suTi = suT[i]; for (int j = 0; j < suTi.length; ++j) { suTi[j] *= a; } } - pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false)); + pseudoInverse = v.multiply(new Array2DRowRealMatrix(suT, false)); this.nonSingular = nonSingular; } - /** Solve the linear equation A × X = B in least square sense. - *The m×n matrix A may not be square, the solution X is - * such that ||A × X - B|| is minimal.
- * @param b right-hand side of the equation A × X = B + /** + * Solve the linear equation A × X = B in least square sense. + *+ * The m×n matrix A may not be square, the solution X is such that + * ||A × X - B|| is minimal. + *
+ * @param b + * right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B - * @exception IllegalArgumentException if matrices dimensions don't match + * @exception IllegalArgumentException + * if matrices dimensions don't match */ - public double[] solve(final double[] b) - throws IllegalArgumentException { + public double[] solve(final double[] b) throws IllegalArgumentException { return pseudoInverse.operate(b); } - /** Solve the linear equation A × X = B in least square sense. - *The m×n matrix A may not be square, the solution X is - * such that ||A × X - B|| is minimal.
- * @param b right-hand side of the equation A × X = B + /** + * Solve the linear equation A × X = B in least square sense. + *+ * The m×n matrix A may not be square, the solution X is such that + * ||A × X - B|| is minimal. + *
+ * @param b + * right-hand side of the equation A × X = B * @return a vector X that minimizes the two norm of A × X - B - * @exception IllegalArgumentException if matrices dimensions don't match + * @exception IllegalArgumentException + * if matrices dimensions don't match */ public RealVector solve(final RealVector b) - throws IllegalArgumentException { + throws IllegalArgumentException { return pseudoInverse.operate(b); } - /** Solve the linear equation A × X = B in least square sense. - *The m×n matrix A may not be square, the solution X is - * such that ||A × X - B|| is minimal.
- * @param b right-hand side of the equation A × X = B + /** + * Solve the linear equation A × X = B in least square sense. + *+ * The m×n matrix A may not be square, the solution X is such that + * ||A × X - B|| is minimal. + *
+ * @param b + * right-hand side of the equation A × X = B * @return a matrix X that minimizes the two norm of A × X - B - * @exception IllegalArgumentException if matrices dimensions don't match + * @exception IllegalArgumentException + * if matrices dimensions don't match */ public RealMatrix solve(final RealMatrix b) - throws IllegalArgumentException { + throws IllegalArgumentException { return pseudoInverse.multiply(b); } @@ -432,7 +363,8 @@ public class SingularValueDecompositionImpl implements SingularValueDecompositio return nonSingular; } - /** Get the pseudo-inverse of the decomposed matrix. + /** + * Get the pseudo-inverse of the decomposed matrix. * @return inverse matrix */ public RealMatrix getInverse() { diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml index d02fbdda1..31c91aac8 100644 --- a/src/site/xdoc/changes.xml +++ b/src/site/xdoc/changes.xml @@ -39,6 +39,12 @@ The