SchurTransformer.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.linear;

  25. import org.hipparchus.exception.LocalizedCoreFormats;
  26. import org.hipparchus.exception.MathIllegalArgumentException;
  27. import org.hipparchus.exception.MathIllegalStateException;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.Precision;

  31. /**
  32.  * Class transforming a general real matrix to Schur form.
  33.  * <p>A m &times; m matrix A can be written as the product of three matrices: A = P
  34.  * &times; T &times; P<sup>T</sup> with P an orthogonal matrix and T an quasi-triangular
  35.  * matrix. Both P and T are m &times; m matrices.</p>
  36.  * <p>Transformation to Schur form is often not a goal by itself, but it is an
  37.  * intermediate step in more general decomposition algorithms like
  38.  * {@link EigenDecompositionSymmetric eigen decomposition}. This class is therefore
  39.  * intended for expert use. As a consequence of this explicitly limited scope,
  40.  * many methods directly returns references to internal arrays, not copies.</p>
  41.  * <p>This class is based on the method hqr2 in class EigenvalueDecomposition
  42.  * from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
  43.  *
  44.  * @see <a href="http://mathworld.wolfram.com/SchurDecomposition.html">Schur Decomposition - MathWorld</a>
  45.  * @see <a href="http://en.wikipedia.org/wiki/Schur_decomposition">Schur Decomposition - Wikipedia</a>
  46.  * @see <a href="http://en.wikipedia.org/wiki/Householder_transformation">Householder Transformations</a>
  47.  */
  48. public class SchurTransformer {
  49.     /** Maximum allowed iterations for convergence of the transformation. */
  50.     private static final int MAX_ITERATIONS = 100;

  52.     /** P matrix. */
  53.     private final double[][] matrixP;
  54.     /** T matrix. */
  55.     private final double[][] matrixT;
  56.     /** Cached value of P. */
  57.     private RealMatrix cachedP;
  58.     /** Cached value of T. */
  59.     private RealMatrix cachedT;
  60.     /** Cached value of PT. */
  61.     private RealMatrix cachedPt;

  63.     /** Epsilon criteria. */
  64.     private final double epsilon;

  66.     /**
  67.      * Build the transformation to Schur form of a general real matrix.
  68.      *
  69.      * @param matrix matrix to transform
  70.      * @throws MathIllegalArgumentException if the matrix is not square
  71.      */
  72.     public SchurTransformer(final RealMatrix matrix) {
  73.         /** Epsilon criteria taken from JAMA code (originally was 2^-52). */
  74.         this(matrix, Precision.EPSILON);
  75.     }

  77.     /**
  78.      * Build the transformation to Schur form of a general real matrix.
  79.      *
  80.      * @param matrix matrix to transform
  81.      * @param epsilon convergence criteria
  82.      * @throws MathIllegalArgumentException if the matrix is not square
  83.      * @since 3.0
  84.      */
  85.     public SchurTransformer(final RealMatrix matrix, final double epsilon) {
  86.         if (!matrix.isSquare()) {
  87.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_SQUARE_MATRIX,
  88.                                                    matrix.getRowDimension(), matrix.getColumnDimension());
  89.         }
  90.         this.epsilon = epsilon;

  92.         HessenbergTransformer transformer = new HessenbergTransformer(matrix);
  93.         matrixT = transformer.getH().getData();
  94.         matrixP = transformer.getP().getData();
  95.         cachedT = null;
  96.         cachedP = null;
  97.         cachedPt = null;

  99.         // transform matrix
  100.         transform();
  101.     }

  103.     /**
  104.      * Returns the matrix P of the transform.
  105.      * <p>P is an orthogonal matrix, i.e. its inverse is also its transpose.</p>
  106.      *
  107.      * @return the P matrix
  108.      */
  109.     public RealMatrix getP() {
  110.         if (cachedP == null) {
  111.             cachedP = MatrixUtils.createRealMatrix(matrixP);
  112.         }
  113.         return cachedP;
  114.     }

  116.     /**
  117.      * Returns the transpose of the matrix P of the transform.
  118.      * <p>P is an orthogonal matrix, i.e. its inverse is also its transpose.</p>
  119.      *
  120.      * @return the transpose of the P matrix
  121.      */
  122.     public RealMatrix getPT() {
  123.         if (cachedPt == null) {
  124.             cachedPt = getP().transpose();
  125.         }

  127.         // return the cached matrix
  128.         return cachedPt;
  129.     }

  131.     /**
  132.      * Returns the quasi-triangular Schur matrix T of the transform.
  133.      *
  134.      * @return the T matrix
  135.      */
  136.     public RealMatrix getT() {
  137.         if (cachedT == null) {
  138.             cachedT = MatrixUtils.createRealMatrix(matrixT);
  139.         }

  141.         // return the cached matrix
  142.         return cachedT;
  143.     }

  145.     /**
  146.      * Transform original matrix to Schur form.
  147.      * @throws MathIllegalStateException if the transformation does not converge
  148.      */
  149.     private void transform() {
  150.         final int n = matrixT.length;

  152.         // compute matrix norm
  153.         final double norm = getNorm();

  155.         // shift information
  156.         final ShiftInfo shift = new ShiftInfo();

  158.         // Outer loop over eigenvalue index
  159.         int iteration = 0;
  160.         int iu = n - 1;
  161.         while (iu >= 0) {

  163.             // Look for single small sub-diagonal element
  164.             final int il = findSmallSubDiagonalElement(iu, norm);

  166.             // Check for convergence
  167.             if (il == iu) {
  168.                 // One root found
  169.                 matrixT[iu][iu] += shift.exShift;
  170.                 iu--;
  171.                 iteration = 0;
  172.             } else if (il == iu - 1) {
  173.                 // Two roots found
  174.                 double p = (matrixT[iu - 1][iu - 1] - matrixT[iu][iu]) / 2.0;
  175.                 double q = p * p + matrixT[iu][iu - 1] * matrixT[iu - 1][iu];
  176.                 matrixT[iu][iu] += shift.exShift;
  177.                 matrixT[iu - 1][iu - 1] += shift.exShift;

  179.                 if (q >= 0) {
  180.                     double z = FastMath.sqrt(FastMath.abs(q));
  181.                     if (p >= 0) {
  182.                         z = p + z;
  183.                     } else {
  184.                         z = p - z;
  185.                     }
  186.                     final double x = matrixT[iu][iu - 1];
  187.                     final double s = FastMath.abs(x) + FastMath.abs(z);
  188.                     p = x / s;
  189.                     q = z / s;
  190.                     final double r = FastMath.sqrt(p * p + q * q);
  191.                     p /= r;
  192.                     q /= r;

  194.                     // Row modification
  195.                     for (int j = iu - 1; j < n; j++) {
  196.                         z = matrixT[iu - 1][j];
  197.                         matrixT[iu - 1][j] = q * z + p * matrixT[iu][j];
  198.                         matrixT[iu][j] = q * matrixT[iu][j] - p * z;
  199.                     }

  201.                     // Column modification
  202.                     for (int i = 0; i <= iu; i++) {
  203.                         z = matrixT[i][iu - 1];
  204.                         matrixT[i][iu - 1] = q * z + p * matrixT[i][iu];
  205.                         matrixT[i][iu] = q * matrixT[i][iu] - p * z;
  206.                     }

  208.                     // Accumulate transformations
  209.                     for (int i = 0; i <= n - 1; i++) {
  210.                         z = matrixP[i][iu - 1];
  211.                         matrixP[i][iu - 1] = q * z + p * matrixP[i][iu];
  212.                         matrixP[i][iu] = q * matrixP[i][iu] - p * z;
  213.                     }
  214.                 }
  215.                 iu -= 2;
  216.                 iteration = 0;
  217.             } else {
  218.                 // No convergence yet
  219.                 computeShift(il, iu, iteration, shift);

  221.                 // stop transformation after too many iterations
  222.                 ++iteration;
  223.                 if (iteration > MAX_ITERATIONS) {
  224.                     throw new MathIllegalStateException(LocalizedCoreFormats.CONVERGENCE_FAILED,
  225.                                                         MAX_ITERATIONS);
  226.                 }

  228.                 // the initial houseHolder vector for the QR step
  229.                 final double[] hVec = new double[3];

  231.                 final int im = initQRStep(il, iu, shift, hVec);
  232.                 performDoubleQRStep(il, im, iu, shift, hVec, norm);
  233.             }
  234.         }
  235.     }

  237.     /**
  238.      * Computes the L1 norm of the (quasi-)triangular matrix T.
  239.      *
  240.      * @return the L1 norm of matrix T
  241.      */
  242.     private double getNorm() {
  243.         double norm = 0.0;
  244.         for (int i = 0; i < matrixT.length; i++) {
  245.             // as matrix T is (quasi-)triangular, also take the sub-diagonal element into account
  246.             for (int j = FastMath.max(i - 1, 0); j < matrixT.length; j++) {
  247.                 norm += FastMath.abs(matrixT[i][j]);
  248.             }
  249.         }
  250.         return norm;
  251.     }

  253.     /**
  254.      * Find the first small sub-diagonal element and returns its index.
  255.      *
  256.      * @param startIdx the starting index for the search
  257.      * @param norm the L1 norm of the matrix
  258.      * @return the index of the first small sub-diagonal element
  259.      */
  260.     private int findSmallSubDiagonalElement(final int startIdx, final double norm) {
  261.         int l = startIdx;
  262.         while (l > 0) {
  263.             double s = FastMath.abs(matrixT[l - 1][l - 1]) + FastMath.abs(matrixT[l][l]);
  264.             if (s == 0.0) {
  265.                 s = norm;
  266.             }
  267.             if (FastMath.abs(matrixT[l][l - 1]) < epsilon * s) {
  268.                 break;
  269.             }
  270.             l--;
  271.         }
  272.         return l;
  273.     }

  275.     /**
  276.      * Compute the shift for the current iteration.
  277.      *
  278.      * @param l the index of the small sub-diagonal element
  279.      * @param idx the current eigenvalue index
  280.      * @param iteration the current iteration
  281.      * @param shift holder for shift information
  282.      */
  283.     private void computeShift(final int l, final int idx, final int iteration, final ShiftInfo shift) {
  284.         // Form shift
  285.         shift.x = matrixT[idx][idx];
  286.         shift.y = shift.w = 0.0;
  287.         if (l < idx) {
  288.             shift.y = matrixT[idx - 1][idx - 1];
  289.             shift.w = matrixT[idx][idx - 1] * matrixT[idx - 1][idx];
  290.         }

  292.         // Wilkinson's original ad hoc shift
  293.         if (iteration == 10) {
  294.             shift.exShift += shift.x;
  295.             for (int i = 0; i <= idx; i++) {
  296.                 matrixT[i][i] -= shift.x;
  297.             }
  298.             final double s = FastMath.abs(matrixT[idx][idx - 1]) + FastMath.abs(matrixT[idx - 1][idx - 2]);
  299.             shift.x = 0.75 * s;
  300.             shift.y = 0.75 * s;
  301.             shift.w = -0.4375 * s * s;
  302.         }

  304.         // MATLAB's new ad hoc shift
  305.         if (iteration == 30) {
  306.             double s = (shift.y - shift.x) / 2.0;
  307.             s = s * s + shift.w;
  308.             if (s > 0.0) {
  309.                 s = FastMath.sqrt(s);
  310.                 if (shift.y < shift.x) {
  311.                     s = -s;
  312.                 }
  313.                 s = shift.x - shift.w / ((shift.y - shift.x) / 2.0 + s);
  314.                 for (int i = 0; i <= idx; i++) {
  315.                     matrixT[i][i] -= s;
  316.                 }
  317.                 shift.exShift += s;
  318.                 shift.x = shift.y = shift.w = 0.964;
  319.             }
  320.         }
  321.     }

  323.     /**
  324.      * Initialize the householder vectors for the QR step.
  325.      *
  326.      * @param il the index of the small sub-diagonal element
  327.      * @param iu the current eigenvalue index
  328.      * @param shift shift information holder
  329.      * @param hVec the initial houseHolder vector
  330.      * @return the start index for the QR step
  331.      */
  332.     private int initQRStep(int il, final int iu, final ShiftInfo shift, double[] hVec) {
  333.         // Look for two consecutive small sub-diagonal elements
  334.         int im = iu - 2;
  335.         while (im >= il) {
  336.             final double z = matrixT[im][im];
  337.             final double r = shift.x - z;
  338.             double s = shift.y - z;
  339.             hVec[0] = (r * s - shift.w) / matrixT[im + 1][im] + matrixT[im][im + 1];
  340.             hVec[1] = matrixT[im + 1][im + 1] - z - r - s;
  341.             hVec[2] = matrixT[im + 2][im + 1];

  343.             if (im == il) {
  344.                 break;
  345.             }

  347.             final double lhs = FastMath.abs(matrixT[im][im - 1]) * (FastMath.abs(hVec[1]) + FastMath.abs(hVec[2]));
  348.             final double rhs = FastMath.abs(hVec[0]) * (FastMath.abs(matrixT[im - 1][im - 1]) +
  349.                                                         FastMath.abs(z) +
  350.                                                         FastMath.abs(matrixT[im + 1][im + 1]));

  352.             if (lhs < epsilon * rhs) {
  353.                 break;
  354.             }
  355.             im--;
  356.         }

  358.         return im;
  359.     }

  361.     /**
  362.      * Perform a double QR step involving rows l:idx and columns m:n
  363.      *
  364.      * @param il the index of the small sub-diagonal element
  365.      * @param im the start index for the QR step
  366.      * @param iu the current eigenvalue index
  367.      * @param shift shift information holder
  368.      * @param hVec the initial houseHolder vector
  369.      * @param norm matrix norm
  370.      */
  371.     private void performDoubleQRStep(final int il, final int im, final int iu,
  372.                                      final ShiftInfo shift, final double[] hVec,
  373.                                      final double norm) {

  375.         final int n = matrixT.length;
  376.         double p = hVec[0];
  377.         double q = hVec[1];
  378.         double r = hVec[2];

  380.         for (int k = im; k <= iu - 1; k++) {
  381.             boolean notlast = k != (iu - 1);
  382.             if (k != im) {
  383.                 p = matrixT[k][k - 1];
  384.                 q = matrixT[k + 1][k - 1];
  385.                 r = notlast ? matrixT[k + 2][k - 1] : 0.0;
  386.                 shift.x = FastMath.abs(p) + FastMath.abs(q) + FastMath.abs(r);
  387.                 if (Precision.equals(shift.x, 0.0, epsilon * norm)) {
  388.                     continue;
  389.                 }
  390.                 p /= shift.x;
  391.                 q /= shift.x;
  392.                 r /= shift.x;
  393.             }
  394.             double s = FastMath.sqrt(p * p + q * q + r * r);
  395.             if (p < 0.0) {
  396.                 s = -s;
  397.             }
  398.             if (s != 0.0) {
  399.                 if (k != im) {
  400.                     matrixT[k][k - 1] = -s * shift.x;
  401.                 } else if (il != im) {
  402.                     matrixT[k][k - 1] = -matrixT[k][k - 1];
  403.                 }
  404.                 p += s;
  405.                 shift.x = p / s;
  406.                 shift.y = q / s;
  407.                 double z = r / s;
  408.                 q /= p;
  409.                 r /= p;

  411.                 // Row modification
  412.                 for (int j = k; j < n; j++) {
  413.                     p = matrixT[k][j] + q * matrixT[k + 1][j];
  414.                     if (notlast) {
  415.                         p += r * matrixT[k + 2][j];
  416.                         matrixT[k + 2][j] -= p * z;
  417.                     }
  418.                     matrixT[k][j] -= p * shift.x;
  419.                     matrixT[k + 1][j] -= p * shift.y;
  420.                 }

  422.                 // Column modification
  423.                 for (int i = 0; i <= FastMath.min(iu, k + 3); i++) {
  424.                     p = shift.x * matrixT[i][k] + shift.y * matrixT[i][k + 1];
  425.                     if (notlast) {
  426.                         p += z * matrixT[i][k + 2];
  427.                         matrixT[i][k + 2] -= p * r;
  428.                     }
  429.                     matrixT[i][k] -= p;
  430.                     matrixT[i][k + 1] -= p * q;
  431.                 }

  433.                 // Accumulate transformations
  434.                 final int high = matrixT.length - 1;
  435.                 for (int i = 0; i <= high; i++) {
  436.                     p = shift.x * matrixP[i][k] + shift.y * matrixP[i][k + 1];
  437.                     if (notlast) {
  438.                         p += z * matrixP[i][k + 2];
  439.                         matrixP[i][k + 2] -= p * r;
  440.                     }
  441.                     matrixP[i][k] -= p;
  442.                     matrixP[i][k + 1] -= p * q;
  443.                 }
  444.             }  // (s != 0)
  445.         }  // k loop

  447.         // clean up pollution due to round-off errors
  448.         for (int i = im + 2; i <= iu; i++) {
  449.             matrixT[i][i-2] = 0.0;
  450.             if (i > im + 2) {
  451.                 matrixT[i][i-3] = 0.0;
  452.             }
  453.         }
  454.     }

  456.     /**
  457.      * Internal data structure holding the current shift information.
  458.      * Contains variable names as present in the original JAMA code.
  459.      */
  460.     private static class ShiftInfo {
  461.         // CHECKSTYLE: stop all

  463.         /** x shift info */
  464.         double x;
  465.         /** y shift info */
  466.         double y;
  467.         /** w shift info */
  468.         double w;
  469.         /** Indicates an exceptional shift. */
  470.         double exShift;

  472.         // CHECKSTYLE: resume all
  473.     }
  474. }
Created with JaCoCo 0.8.12.202403310830