Class RealSchur

java.lang.Object
org.flag4j.linalg.decompositions.Decomposition<Matrix>
org.flag4j.linalg.decompositions.schur.Schur<Matrix,double[]>
org.flag4j.linalg.decompositions.schur.RealSchur
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.
public class RealSchur extends Schur<Matrix,double[]>

Instanced of this class can be used for computing the Schur decomposition of a real dense square matrix.

The Schur decomposition decomposes a given square matrix \( A \) into: \[ A = UTU^{T} \] where \( U \) is an orthogonal 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 an instance of RealSchur.
  2. Call decompose(Matrix) to perform the factorization.
  3. Retrieve the resulting matrices using Schur.getU() and Schur.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 double
    currentFactor
    Stores the scalar factor \( \alpha \) for use in computation of the Householder reflector \( P = I - \alpha vv^{T} \).
    protected double
    norm
    For computing the norm of a column for use when computing Householder reflectors.

    Fields inherited from class org.flag4j.linalg.decompositions.schur.Schur

    balancer, checkFinite, computeU, DEFAULT_EXCEPTIONAL_ITERS, DEFAULT_MAX_ITERS_FACTOR, exceptionalThreshold, hess, householderVector, iHigh, iLow, maxIterations, maxIterationsFactor, numExceptional, numRows, rng, shiftCol, sinceLastExceptional, T, temp, U, workArray

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

    hasDecomposed

  • Constructor Summary

    Constructors
    Constructor
    Description
    RealSchur()
    Creates a decomposer to compute the Schur decomposition for a real dense matrix.
    RealSchur(boolean computeU)
    Creates a decomposer to compute the Schur decomposition for a real dense matrix where the \( U \) matrix may or may not be computed.
    RealSchur(boolean computeU, long seed)
    Creates a decomposer to compute the Schur decomposition for a real dense matrix.
    RealSchur(long seed)
    Creates a decomposer to compute the Schur decomposition for a real dense matrix.

  • Method Summary

    Modifier and Type
    Method
    Description
    protected void
    applyDoubleShiftReflector(int i, boolean set)
    Applies reflector for the double shift.
    protected void
    applyReflector(int i, int shiftSize)
    Applies the constructed Householder reflector which has been constructed for the given shift size.
    protected void
    applySingleShiftReflector(int i, boolean set)
    Applies reflector for the double shift.
    protected 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 void
    checkFinite(Matrix src)
    Ensures that src only contains finite values.
    protected double
    computeExceptionalShift(int k)
    Computes a random shift to help the QR algorithm converge if it gets stuck.
    protected void
    computeImplicitDoubleShift(int workEnd)
    Computes the shifts for a Francis double shift iteration.
    protected void
    computeImplicitSingleShift(int k, double shift)
    Computes the non-zero data of the first column for the single shifted QR algorithm.
    RealSchur
    decompose(Matrix src)
    Computes the Schur decomposition of the input matrix.
    RealSchur
    enforceFinite(boolean enforceFinite)
    Sets flag indicating if a check should be made to ensure the matrix being decomposed only contains finite values.
    protected boolean
    makeReflector(int i, double p1, double p2)
    Constructs a householder reflector given specified values for a column to apply the reflector to.
    protected boolean
    makeReflector(int i, double p1, double p2, double p3)
    Constructs a householder reflector given specified values for a column to apply the reflector to.
    protected void
    performDoubleShift(int workEnd)
    Performs a full iteration of the Francis implicit double shifted QR algorithm (this includes the bulge chase).
    protected 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.
    protected void
    performSingleShift(int workEnd, double shift)
    Performs a full iteration of the implicit single shifted QR algorithm (this includes the bulge chase).
    CMatrix[]
    real2ComplexSchur()
    Converts the real schur form computed in the last decomposition to the complex Schur form.
    RealSchur
    setExceptionalThreshold(int exceptionalThreshold)
    Sets the number of iterations of the QR algorithm to perform without deflation before performing a random shift.
    RealSchur
    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
    setUpArrays()
    Initializes temporary work arrays to be used in the decomposition.
    protected 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.schur.Schur

    decomposeBase, getT, getU, setUp

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

    ensureHasDecomposed

    Methods inherited from class java.lang.Object

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

  • Field Details

    • norm

      protected double norm
      For computing the norm of a column for use when computing Householder reflectors.

    • currentFactor

      protected double currentFactor
      Stores the scalar factor \( \alpha \) for use in computation of the Householder reflector \( P = I - \alpha vv^{T} \).

  • Constructor Details

    • RealSchur

      public RealSchur()

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

      Note: This decomposer may use random numbers during the decomposition. If reproducible results are needed, set the seed for the pseudo-random number generator using RealSchur(long)

    • RealSchur

      public RealSchur(boolean computeU)

      Creates a decomposer to compute the Schur decomposition for a real dense matrix where the \( U \) matrix may or may not be computed.

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

      By default, if a constructor with no computeU parameter is called, \( U \) will be computed.

      Note: This decomposer may use random numbers during the decomposition. If reproducible results are needed, set the seed for the pseudo-random number generator using RealSchur(boolean, long)

      Parameters:
      computeU - Flag indicating if the unitary \( U \) matrix should be computed for the Schur decomposition. If true, \( U \) will be computed. If false, \( U \) will not be computed.

    • RealSchur

      public RealSchur(long seed)
      Creates a decomposer to compute the Schur decomposition for a real dense matrix.
      Parameters:
      seed - Seed to use for pseudo-random number generator when computing exceptional shifts during the QR algorithm.

    • RealSchur

      public RealSchur(boolean computeU, long seed)
      Creates a decomposer to compute the Schur decomposition for a real dense matrix.
      Parameters:
      seed - Seed to use for pseudo-random number generator when computing exceptional shifts during the QR algorithm.

  • Method Details

    • setExceptionalThreshold

      public RealSchur setExceptionalThreshold(int exceptionalThreshold)
      Description copied from class: Schur

      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 Schur.DEFAULT_EXCEPTIONAL_ITERS

      Overrides:
      setExceptionalThreshold in class Schur<Matrix,double[]>
      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.

    • setMaxIterationFactor

      public RealSchur setMaxIterationFactor(int maxIterationFactor)
      Description copied from class: Schur

      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 Schur.DEFAULT_MAX_ITERS_FACTOR).

      Overrides:
      setMaxIterationFactor in class Schur<Matrix,double[]>
      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.

    • enforceFinite

      public RealSchur enforceFinite(boolean enforceFinite)
      Description copied from class: Schur

      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.

      Overrides:
      enforceFinite in class Schur<Matrix,double[]>
      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.

    • unbalance

      protected 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.

      Specified by:
      unbalance in class Schur<Matrix,double[]>

    • decompose

      public RealSchur decompose(Matrix src)

      Computes the Schur decomposition of the input matrix.

      Specified by:
      decompose in class Decomposition<Matrix>
      Parameters:
      src - The source matrix to decompose.
      Returns:
      A reference to this decomposer.

    • setUpArrays

      protected void setUpArrays()
      Initializes temporary work arrays to be used in the decomposition.
      Specified by:
      setUpArrays in class Schur<Matrix,double[]>

    • performExceptionalShift

      protected 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.
      Specified by:
      performExceptionalShift in class Schur<Matrix,double[]>
      Parameters:
      workEnd - The ending row (inclusive) of the current active working block.

    • computeExceptionalShift

      protected double computeExceptionalShift(int k)
      Computes a random shift to help the QR algorithm converge if it gets stuck.
      Parameters:
      k - The current size of the working matrix. Specifically, the index of the lower right value in the working matrix is (k, k).
      Returns:
      A shift in a random direction which has the same magnitude as the elements in the matrix.

    • computeImplicitSingleShift

      protected void computeImplicitSingleShift(int k, double shift)
      Computes the non-zero data of the first column for the single shifted QR algorithm.
      Parameters:
      k - Size of current working matrix.
      shift - The shift to use.

    • performSingleShift

      protected void performSingleShift(int workEnd, double shift)
      Performs a full iteration of the implicit single shifted QR algorithm (this includes the bulge chase).
      Parameters:
      workEnd - The ending row (inclusive) of the current active working block.
      shift - The shift to use in the implicit single shifted QR algorithm.

    • applySingleShiftReflector

      protected void applySingleShiftReflector(int i, boolean set)
      Applies reflector for the double shift. This method can be used to apply either be the reflector constructed for the first column of the shifted matrix, or a reflector being used in the bulge chase of size 2 which arises from the first case.
      Parameters:
      i - The starting row the reflector is being applied to.

    • performDoubleShift

      protected void performDoubleShift(int workEnd)
      Performs a full iteration of the Francis implicit double shifted QR algorithm (this includes the bulge chase).
      Specified by:
      performDoubleShift in class Schur<Matrix,double[]>
      Parameters:
      workEnd - The ending row (inclusive) of the current active working block.

    • computeImplicitDoubleShift

      protected void computeImplicitDoubleShift(int workEnd)
      Computes the shifts for a Francis double shift iteration. Specifically, the shifts are the generalized Rayleigh quotients of degree two.
      Parameters:
      workEnd - The ending row (inclusive) of the current active working block.

    • applyDoubleShiftReflector

      protected void applyDoubleShiftReflector(int i, boolean set)
      Applies reflector for the double shift. This method can be used to apply either be the reflector constructed for the first column of the shifted matrix, or a reflector being used in the bulge chase of size 2 which arises from the first case.
      Parameters:
      i - The starting row the reflector is being applied to.

    • applyReflector

      protected void applyReflector(int i, int shiftSize)
      Applies the constructed Householder reflector which has been constructed for the given shift size.
      Parameters:
      i - The stating row the reflector is being applied to.
      shiftSize - The size of the shift the reflector was constructed for.

    • makeReflector

      protected boolean makeReflector(int i, double p1, double p2, double p3)
      Constructs a householder reflector given specified values for a column to apply the reflector to. This reflector is stored in indices i, i+1, and i+2 of Schur.householderVector.
      Parameters:
      i - Row of working matrix to construct reflector for.
      p1 - First entry to in column to apply reflector to.
      p2 - Second entry in column to apply reflector to.
      p3 - Third entry in column to apply reflector to.
      Returns:
      True if a reflector needs to be constructed to return matrix to upper Hessenburg form. False if column is already in the correct form.

    • makeReflector

      protected boolean makeReflector(int i, double p1, double p2)
      Constructs a householder reflector given specified values for a column to apply the reflector to. This reflector is stored in indices i and i+1 of Schur.householderVector.
      Parameters:
      i - Row of working matrix to construct reflector for.
      p1 - First entry to in column to apply reflector to.
      p2 - Second entry in column to apply reflector to.
      Returns:
      True if a reflector needs to be constructed to return matrix to upper Hessenburg form. False if column is already in the correct form.

    • checkConvergence

      protected 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.
      Specified by:
      checkConvergence in class Schur<Matrix,double[]>
      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 void checkFinite(Matrix src)
      Ensures that src only contains finite values.
      Specified by:
      checkFinite in class Schur<Matrix,double[]>
      Parameters:
      src - Matrix of interest.
      Throws:
      IllegalArgumentException - If src does not contain only finite values.

    • real2ComplexSchur

      public CMatrix[] real2ComplexSchur()

      Converts the real schur form computed in the last decomposition to the complex Schur form.

      That is, converts the real block upper triangular Schur matrix to a complex valued properly upper triangular matrix. If the unitary transformation matrix \( U \) was computed, the transformations will also be updated accordingly.

      This method was adapted from the code given by scipy.linalg.rsf2csf (v1.12.0).

      Returns:
      An array of length 2 containing the complex Schur matrix \( T \) from the last decomposition, and if computed, the complex unitary transformation matrix \( U \) from the decomposition. If \( U \) was not computed, then the arrays second value will be null.