Class Schur<T extends MatrixMixin<T,?,?,?>,U>

java.lang.Object
org.flag4j.linalg.decompositions.Decomposition<T>
org.flag4j.linalg.decompositions.schur.Schur<T,U>
Type Parameters:
T - The type of matrix to be decomposed.
U - The type for the internal storage data structure of the matrix to be decomposed.
Direct Known Subclasses:
ComplexSchur, RealSchur
public abstract class Schur<T extends MatrixMixin<T,?,?,?>,U> extends Decomposition<T>

An abstract base class for computing the Schur decomposition of a square matrix.

The Schur decomposition decomposes a given square matrix \( A \) into: \[ A = UTU^{H} \] where \( U \) is a unitary (or orthogonal for real matrices) matrix \( T \) is a quasi-upper triangular matrix known as the Schur form of \( A \). This means \( T \) is upper triangular except for possibly \( 2\times2 \) blocks along its diagonal, which correspond to complex conjugate pairs of eigenvalues.

The Schur decomposition proceeds by an iterative algorithm with possible random behavior. For reproducibility, constructors support specifying a seed for the pseudo-random number generator.

Usage:

The decomposition workflow typically follows these steps:
  1. Instantiate a con concrete instance of Schur.
  2. Call Decomposition.decompose(MatrixMixin) to perform the factorization.
  3. Retrieve the resulting matrices using getU() and getT().

Efficiency Considerations:

If eigenvectors are not required, setting computeU = false may improve performance.

This class was inspired by code from the EJML library and the description of the Francis implicit double shifted QR algorithm from Fundamentals of Matrix Computations 3rd Edition by David S. Watkins.

See Also:

  • Field Summary

    Fields
    Modifier and Type
    Field
    Description
    protected Balancer<T>
    balancer
    Balancer to apply a similarity transform to the matrix before the QR-algorithm is executed.
    protected boolean
    checkFinite
    Flag indicating if a check should be made during the decomposition that the working matrix contains only finite values.
    protected final boolean
    computeU
    Flag indicating if the orthogonal matrix \( U \) in the Schur decomposition should be computed.
    protected final int
    DEFAULT_EXCEPTIONAL_ITERS
    Default number of iterations to apply before doing an exceptional shift.
    protected final int
    DEFAULT_MAX_ITERS_FACTOR
    Default factor for computing the maximum number of iterations to perform.
    protected int
    exceptionalThreshold
    Number of iterations to run without deflation before an exceptional shift is done.
    protected UnitaryDecomposition<T,U>
    hess
    Decomposer to compute the Hessenburg decomposition as a setup step for the implicit double step QR algorithm.
    protected U
    householderVector
    Stores the vector \( v \) in the Householder reflector \( P = I - \alpha vv^{T} \).
    protected int
    iHigh
    The upper bound (exclusive) of the row/column indices of the block to reduce to Schur form.
    protected int
    iLow
    The lower bound (inclusive) of the row/column indices of the block to reduce to Schur form.
    protected int
    maxIterations
    Maximum number of iterations to run QR algorithm for.
    protected int
    maxIterationsFactor
    Factor for computing the maximum number of iterations to run the QR algorithm for.
    protected int
    numExceptional
    The total number of exceptional shifts that have been used during the decomposition.
    protected int
    numRows
    Stores the number of rows in the matrix being decomposed (after balancing).
    protected final RandomComplex
    rng
    Random number generator to be used when computing a random exceptional shift.
    protected U
    shiftCol
    Stores the non-zero data of the first column of the shifted matrix \( \left(A- \rho _{1}I\right)\left(A-\rho _{2} I\right) \) where \( \rho _{1} \) and \( \rho _{2} \) are the two shifts.
    protected int
    sinceLastExceptional
    The number of iterations run in the QR algorithm without deflating or performing an exceptional shift.
    protected T
    T
    For storing the (possibly block) upper triangular matrix \( T \) in the Schur decomposition.
    protected U
    temp
    Array for holding temporary values when computing the shifts.
    protected T
    U
    For storing the unitary \( U \) matrix in the Schur decomposition.
    protected U
    workArray
    An array for storing temporary values along the colum of a matrix when applying Householder reflectors.

    Fields inherited from class org.flag4j.linalg.decompositions.Decomposition

    hasDecomposed

  • Constructor Summary

    Constructors
    Modifier
    Constructor
    Description
    protected
    Schur(boolean computeU, RandomComplex rng, UnitaryDecomposition<T,U> hess, Balancer<T> balancer)
    Creates a decomposer to compute the Schur decomposition for a real dense matrix.

  • Method Summary

    Modifier and Type
    Method
    Description
    protected abstract int
    checkConvergence(int workEnd)
    Checks for convergence of lower \( 2\times2 \) sub-matrix within working matrix to upper triangular or block upper triangular form.
    protected abstract void
    checkFinite(T src)
    Ensures that src only contains finite values.
    protected void
    decomposeBase(T src)
    Computes the Schur decomposition of the input matrix.
    Schur<T,U>
    enforceFinite(boolean enforceFinite)
    Sets flag indicating if a check should be made to ensure the matrix being decomposed only contains finite values.
    T
    getT()
    Gets the upper, or possibly block-upper, triangular Schur matrix \( T \) from the Schur decomposition
    T
    getU()
    Gets the unitary matrix \( U \) from the Schur decomposition containing the Schur vectors as its columns.
    protected abstract void
    performDoubleShift(int workEnd)
    Performs a full iteration of the Francis implicit double shifted QR algorithm (this includes the bulge chase).
    protected abstract void
    performExceptionalShift(int workEnd)
    Performs a full iteration of the single shifted QR algorithm (this includes the bulge chase) where the shift is chosen to be a random value with the same magnitude as the lower right element of the working matrix.
    Schur<T,U>
    setExceptionalThreshold(int exceptionalThreshold)
    Sets the number of iterations of the QR algorithm to perform without deflation before performing a random shift.
    Schur<T,U>
    setMaxIterationFactor(int maxIterationFactor)
    Specify maximum iteration factor for computing the total number of iterations to run the QR algorithm for when computing the decomposition.
    protected void
    setUp(T src)
    Performs basic setup and initializes data structures to be used in the decomposition.
    protected abstract void
    setUpArrays()
    Initializes temporary work arrays to be used in the decomposition.
    protected abstract void
    unbalance()
    Reverts the scaling and permutations applied during the balancing step to obtain the correct form.

    Methods inherited from class org.flag4j.linalg.decompositions.Decomposition

    decompose, ensureHasDecomposed

    Methods inherited from class java.lang.Object

    clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

  • Field Details

    • rng

      protected final RandomComplex rng
      Random number generator to be used when computing a random exceptional shift.

    • DEFAULT_EXCEPTIONAL_ITERS

      protected final int DEFAULT_EXCEPTIONAL_ITERS
      Default number of iterations to apply before doing an exceptional shift.
      See Also:

    • DEFAULT_MAX_ITERS_FACTOR

      protected final int DEFAULT_MAX_ITERS_FACTOR
      Default factor for computing the maximum number of iterations to perform.
      See Also:

    • T

      protected T extends MatrixMixin<T,?,?,?> T
      For storing the (possibly block) upper triangular matrix \( T \) in the Schur decomposition.

    • U

      protected T extends MatrixMixin<T,?,?,?> U
      For storing the unitary \( U \) matrix in the Schur decomposition.

    • hess

      protected UnitaryDecomposition<T extends MatrixMixin<T,?,?,?>,U> hess
      Decomposer to compute the Hessenburg decomposition as a setup step for the implicit double step QR algorithm.

    • balancer

      protected Balancer<T extends MatrixMixin<T,?,?,?>> balancer

      Balancer to apply a similarity transform to the matrix before the QR-algorithm is executed. This similarity transform consists of permuting rows/columns to isolate decoupled eigenvalues then scaling the rows and columns of the matrix

      This is done to attempt to improve the conditioning of the eigen-problem.

    • iLow

      protected int iLow
      The lower bound (inclusive) of the row/column indices of the block to reduce to Schur form.

    • iHigh

      protected int iHigh
      The upper bound (exclusive) of the row/column indices of the block to reduce to Schur form.

    • numRows

      protected int numRows
      Stores the number of rows in the matrix being decomposed (after balancing).

    • householderVector

      protected U householderVector
      Stores the vector \( v \) in the Householder reflector \( P = I - \alpha vv^{T} \).

    • shiftCol

      protected U shiftCol
      Stores the non-zero data of the first column of the shifted matrix \( \left(A- \rho _{1}I\right)\left(A-\rho _{2} I\right) \) where \( \rho _{1} \) and \( \rho _{2} \) are the two shifts.

    • workArray

      protected U workArray
      An array for storing temporary values along the colum of a matrix when applying Householder reflectors. This can help improve cache performance when applying the reflector.

    • temp

      protected U temp
      Array for holding temporary values when computing the shifts.

    • maxIterationsFactor

      protected int maxIterationsFactor
      Factor for computing the maximum number of iterations to run the QR algorithm for.

    • maxIterations

      protected int maxIterations
      Maximum number of iterations to run QR algorithm for.

    • exceptionalThreshold

      protected int exceptionalThreshold
      Number of iterations to run without deflation before an exceptional shift is done.

    • sinceLastExceptional

      protected int sinceLastExceptional
      The number of iterations run in the QR algorithm without deflating or performing an exceptional shift.

    • numExceptional

      protected int numExceptional
      The total number of exceptional shifts that have been used during the decomposition.

    • checkFinite

      protected boolean checkFinite
      Flag indicating if a check should be made during the decomposition that the working matrix contains only finite values. If true, an explicit check will be made and an exception will be thrown if non-finite values are found. If false, no check will be made and the floating point arithmetic will carry on with infinities, negative-infinities, and NaNs present.

    • computeU

      protected final boolean computeU
      Flag indicating if the orthogonal matrix \( U \) in the Schur decomposition should be computed.
      • If true, \( U \) will be computed.
      • If false, \( U \) will not be computed. This may improve performance if \( U \) is not required.

  • Constructor Details

    • Schur

      protected Schur(boolean computeU, RandomComplex rng, UnitaryDecomposition<T,U> hess, Balancer<T> balancer)

      Creates a decomposer to compute the Schur decomposition for a real dense matrix.

      If the \( U \) matrix is not needed, passing computeU = false may provide a performance improvement.

      Parameters:
      computeU - Flag indicating if the orthogonal matrix \( U \) in the Schur decomposition should be computed.
      • If true, \( U \) will be computed.
      • If false, \( U \) will not be computed. This may improve performance if \( U \) is not required.
      rng - Random number generator to use when performing random exceptional shifts.
      hess - Decomposer to compute the Hessenburg decomposition as a setup step for the QR algorithm.
      balancer - Balancer which balances the matrix as a preprocessing step to improve the conditioning of the eigenvalue problem.

  • Method Details

    • setExceptionalThreshold

      public Schur<T,U> setExceptionalThreshold(int exceptionalThreshold)

      Sets the number of iterations of the QR algorithm to perform without deflation before performing a random shift.

      That is, if exceptionalThreshold = 10, then at most 10 iterations QR algorithm iterations will be performed. If, by the 10th iteration, no convergence has been detected which allows for deflation, then a QR algorithm iteration will be performed with a random (i.e. exceptional) shift.

      By default, the threshold is set to DEFAULT_EXCEPTIONAL_ITERS

      Parameters:
      exceptionalThreshold - The new exceptional shift threshold. i.e. the number of iterations to perform without deflation before performing an iteration with random shifts.
      Returns:
      A reference to this Schur decomposer.
      Throws:
      IllegalArgumentException - If exceptionalThreshold is not positive.

    • setMaxIterationFactor

      public Schur<T,U> setMaxIterationFactor(int maxIterationFactor)

      Specify maximum iteration factor for computing the total number of iterations to run the QR algorithm for when computing the decomposition. The maximum number of iterations is computed as maxIteration = maxIterationFactor * src.numRows; If the algorithm does not converge within this limit, an exception will be thrown.

      By default, this is computed as maxIterations = DEFAULT_MAX_ITERS_FACTOR * src.numRows; where src is the matrix being decomposed (see DEFAULT_MAX_ITERS_FACTOR).

      Parameters:
      maxIterationFactor - maximum iteration factor for use in computing the total maximum number of iterations to run the QR algorithm for.
      Returns:
      A reference to this Schur decomposer.
      Throws:
      IllegalArgumentException - If maxIterationFactor is not positive.

    • enforceFinite

      public Schur<T,U> enforceFinite(boolean enforceFinite)

      Sets flag indicating if a check should be made to ensure the matrix being decomposed only contains finite values.

      By default, this will be false.

      Parameters:
      enforceFinite - Flag indicating if a check should be made to ensure matrices decomposed by this instance only contain finite values.
      • If true, an explicit check will be made.
      • If false, an explicit check will not be made.
      Returns:
      A reference to this Schur decomposer.

    • getT

      public T getT()
      Gets the upper, or possibly block-upper, triangular Schur matrix \( T \) from the Schur decomposition
      Returns:
      The \( T \) matrix from the Schur decomposition \( A=UTU^{H} \).

    • getU

      public T getU()
      Gets the unitary matrix \( U \) from the Schur decomposition containing the Schur vectors as its columns.
      Returns:
      The \( U \) matrix from the Schur decomposition \( A=UTU^{H} \).

    • decomposeBase

      protected void decomposeBase(T src)

      Computes the Schur decomposition of the input matrix.

      Parameters:
      src - The source matrix to decompose.
      Throws:
      LinearAlgebraException - If the decomposition does not converge within the specified number of max iterations.

    • setUp

      protected void setUp(T src)
      Performs basic setup and initializes data structures to be used in the decomposition.
      Parameters:
      src - The matrix to be decomposed.

    • unbalance

      protected abstract void unbalance()

      Reverts the scaling and permutations applied during the balancing step to obtain the correct form.

      Specifically, this method computes \[ \begin{align*} U :&= PDU \\ &= TU \end{align*} \] where \( P \) and \( D \) are the permutation and scaling matrices respectively from balancing.

    • setUpArrays

      protected abstract void setUpArrays()
      Initializes temporary work arrays to be used in the decomposition.

    • performExceptionalShift

      protected abstract void performExceptionalShift(int workEnd)
      Performs a full iteration of the single shifted QR algorithm (this includes the bulge chase) where the shift is chosen to be a random value with the same magnitude as the lower right element of the working matrix. This can help the QR converge for certain pathological cases where the double shift algorithm oscillates or fails to converge for repeated eigenvalues.
      Parameters:
      workEnd - The ending row (inclusive) of the current active working block.

    • performDoubleShift

      protected abstract void performDoubleShift(int workEnd)
      Performs a full iteration of the Francis implicit double shifted QR algorithm (this includes the bulge chase).
      Parameters:
      workEnd - The ending row (inclusive) of the current active working block.

    • checkConvergence

      protected abstract int checkConvergence(int workEnd)
      Checks for convergence of lower \( 2\times2 \) sub-matrix within working matrix to upper triangular or block upper triangular form. If convergence is found, this will also zero out the values which have converged to near zero.
      Parameters:
      workEnd - The ending row (inclusive) of the current active working block.
      Returns:
      Returns the amount the working matrix size should be deflated. Will be zero if no convergence is detected, one if convergence to upper triangular form is detected and two if convergence to block upper triangular form is detected.

    • checkFinite

      protected abstract void checkFinite(T src)
      Ensures that src only contains finite values.
      Parameters:
      src - Matrix of interest.
      Throws:
      IllegalArgumentException - If src does not contain only finite values.