## org.apache.commons.math.linear
Class EigenDecompositionImpl

    java.lang.Object
      org.apache.commons.math.linear.EigenDecompositionImpl

All Implemented Interfaces:
  EigenDecomposition

    public class EigenDecompositionImpl

extends java.lang.Object
implements EigenDecomposition

Calculates the eigen decomposition of a symmetric matrix.

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 `getD()` is always diagonal and the imaginary values returned `getImagEigenvalue(int)` and `getImagEigenvalues()` are always null.

When called with a `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 `getEigenvector(int)`, `getV()`, `getVT()`, `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).

Since:
  2.0
Version:
  $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
Author:
  Beresford Parlett, University of California, Berkeley, USA (fortran version), Jim Demmel, University of California, Berkeley, USA (fortran version), Inderjit Dhillon, University of Texas, Austin, USA(fortran version), Osni Marques, LBNL/NERSC, USA (fortran version), Christof Voemel, University of California, Berkeley, USA(fortran version)

Nested Class Summary
`private static class` `EigenDecompositionImpl.Solver` Specialized solver.

Field Summary
`private RealMatrix` `cachedD` Cached value of D.
`private RealMatrix` `cachedV` Cached value of V.
`private RealMatrix` `cachedVt` Cached value of Vt.
`private double` `dMin` Minimal value on current state of the diagonal.
`private double` `dMin1` Minimal value on current state of the diagonal, excluding last element.
`private double` `dMin2` Minimal value on current state of the diagonal, excluding last two elements.
`private double` `dN` Last value on current state of the diagonal.
`private double` `dN1` Last but one value on current state of the diagonal.
`private double` `dN2` Last but two on current state of the diagonal.
`private ArrayRealVector[]` `eigenvectors` Eigenvectors.
`private double` `eMin` Min value of off-diagonal elements in current segment.
`private double` `g` Shift ratio with respect to dMin used when tType == 6.
`private double[]` `imagEigenvalues` Imaginary part of the realEigenvalues.
`private double` `lowerSpectra` Lower bound of spectra.
`private double[]` `main` Main diagonal of the tridiagonal matrix.
`private double` `minPivot` Minimum pivot in the Sturm sequence.
`private int` `pingPong` Shift within qd array for ping-pong implementation.
`private double` `qMax` Max value of diagonal elements in current segment.
`private double[]` `realEigenvalues` Real part of the realEigenvalues.
`private double[]` `secondary` Secondary diagonal of the tridiagonal matrix.
`private double` `sigma` Current shift.
`private double` `sigmaLow` Low part of the current shift.
`private double` `splitTolerance` Split tolerance.
`private double[]` `squaredSecondary` Squared secondary diagonal of the tridiagonal matrix.
`private double` `tau` Shift increment to apply.
`private static double` `TOLERANCE` Tolerance.
`private static double` `TOLERANCE_2` Squared tolerance.
`private TriDiagonalTransformer` `transformer` Transformer to tridiagonal (may be null if matrix is already tridiagonal).
`private int` `tType` Type of the last dqds shift.
`private double` `upperSpectra` Upper bound of spectra.
`private double[]` `work` Work array for all decomposition algorithms.

Constructor Summary
`EigenDecompositionImpl(double[] main, double[] secondary, double splitTolerance)` Calculates the eigen decomposition of the given tridiagonal symmetric matrix.
`EigenDecompositionImpl(RealMatrix matrix, double splitTolerance)` Calculates the eigen decomposition of the given symmetric matrix.

Method Summary
`private void` `computeGershgorinCircles()` Compute the Gershgorin circles for all rows.
`private void` `computeShiftIncrement(int start, int end, int deflated)` Compute the shift increment as an estimate of the smallest eigenvalue.
`private java.util.List<java.lang.Integer>` `computeSplits()` Compute splitting points.
`private int` `countEigenValues(double t, int index, int n)` Count the number of realEigenvalues below a point.
`private void` `decompose()` Decompose a tridiagonal symmetric matrix.
`private void` `dqd(int start, int end)` Perform a dqd step.
`private void` `dqds(int start, int end)` Perform a dqds step, using current shift increment.
`private double[]` `eigenvaluesRange(int index, int n)` Compute an interval containing all realEigenvalues of a block.
`private void` `findEigenvalues()` Find the realEigenvalues.
`private ArrayRealVector` `findEigenvector(double eigenvalue, double[] d, double[] l)` Find an eigenvector corresponding to an eigenvalue, using bidiagonals.
`private void` `findEigenVectors()` Find eigenvectors.
`private boolean` `flipIfWarranted(int n, int step)` Flip qd array if warranted.
`RealMatrix` `getD()` Returns the block diagonal matrix D of the decomposition.
`double` `getDeterminant()` Return the determinant of the matrix
`RealVector` `getEigenvector(int i)` Returns a copy of the ith eigenvector of the original matrix.
`double` `getImagEigenvalue(int i)` Returns the imaginary part of the ith eigenvalue of the original matrix.
`double[]` `getImagEigenvalues()` Returns a copy of the imaginary parts of the eigenvalues of the original matrix.
`double` `getRealEigenvalue(int i)` Returns the real part of the ith eigenvalue of the original matrix.
`double[]` `getRealEigenvalues()` Returns a copy of the real parts of the eigenvalues of the original matrix.
`DecompositionSolver` `getSolver()` Get a solver for finding the A × X = B solution in exact linear sense.
`RealMatrix` `getV()` Returns the matrix V of the decomposition.
`RealMatrix` `getVT()` Returns the transpose of the matrix V of the decomposition.
`private int` `goodStep(int start, int end)` Perform one "good" dqd/dqds step.
`private void` `initialSplits(int n)` Perform two iterations with Li's tests for initial splits.
`private boolean` `isSymmetric(RealMatrix matrix)` Check if a matrix is symmetric.
`private void` `ldlTDecomposition(double lambda, int index, int n)` Decompose the shifted tridiagonal matrix T-λI as LDLT.
`private void` `process1RowBlock(int index)` Find eigenvalue in a block with 1 row.
`private void` `process2RowsBlock(int index)` Find realEigenvalues in a block with 2 rows.
`private void` `process3RowsBlock(int index)` Find realEigenvalues in a block with 3 rows.
`private void` `processGeneralBlock(int n)` Find realEigenvalues using dqd/dqds algorithms.
`private void` `progressiveQuotientDifferenceWithShift(double[] d, double[] l, double lambda)` Decompose matrix LDLT - λ I as U-D-U-T.
`private void` `stationaryQuotientDifferenceWithShift(double[] d, double[] l, double lambda)` Decompose matrix LDLT - λ I as L+D+L+T.
`private void` `transformToTridiagonal(RealMatrix matrix)` Transform matrix to tridiagonal.
`private void` `updateSigma(double tau)` Update sigma.

Methods inherited from class java.lang.Object
`clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait`

Field Detail

### TOLERANCE

    private static final double TOLERANCE

  Tolerance.

  See Also:
    Constant Field Values

### TOLERANCE_2

    private static final double TOLERANCE_2

  Squared tolerance.

  See Also:
    Constant Field Values

### splitTolerance

    private double splitTolerance

  Split tolerance.

### main

    private double[] main

  Main diagonal of the tridiagonal matrix.

### secondary

    private double[] secondary

  Secondary diagonal of the tridiagonal matrix.

### squaredSecondary

    private double[] squaredSecondary

  Squared secondary diagonal of the tridiagonal matrix.

### transformer

    private TriDiagonalTransformer transformer

  Transformer to tridiagonal (may be null if matrix is already tridiagonal).

### lowerSpectra

    private double lowerSpectra

  Lower bound of spectra.

### upperSpectra

    private double upperSpectra

  Upper bound of spectra.

### minPivot

    private double minPivot

  Minimum pivot in the Sturm sequence.

### sigma

    private double sigma

  Current shift.

### sigmaLow

    private double sigmaLow

  Low part of the current shift.

### tau

    private double tau

  Shift increment to apply.

### work

    private double[] work

  Work array for all decomposition algorithms.

### pingPong

    private int pingPong

  Shift within qd array for ping-pong implementation.

### qMax

    private double qMax

  Max value of diagonal elements in current segment.

### eMin

    private double eMin

  Min value of off-diagonal elements in current segment.

### tType

    private int tType

  Type of the last dqds shift.

### dMin

    private double dMin

  Minimal value on current state of the diagonal.

### dMin1

    private double dMin1

  Minimal value on current state of the diagonal, excluding last element.

### dMin2

    private double dMin2

  Minimal value on current state of the diagonal, excluding last two elements.

### dN

    private double dN

  Last value on current state of the diagonal.

### dN1

    private double dN1

  Last but one value on current state of the diagonal.

### dN2

    private double dN2

  Last but two on current state of the diagonal.

### g

    private double g

  Shift ratio with respect to dMin used when tType == 6.

### realEigenvalues

    private double[] realEigenvalues

  Real part of the realEigenvalues.

### imagEigenvalues

    private double[] imagEigenvalues

  Imaginary part of the realEigenvalues.

### eigenvectors

    private ArrayRealVector[] eigenvectors

  Eigenvectors.

### cachedV

    private RealMatrix cachedV

  Cached value of V.

### cachedD

    private RealMatrix cachedD

  Cached value of D.

### cachedVt

    private RealMatrix cachedVt

  Cached value of Vt.

Constructor Detail

### EigenDecompositionImpl

    public EigenDecompositionImpl(RealMatrix matrix,
                                  double splitTolerance)
                           throws InvalidMatrixException

  Calculates the eigen decomposition of the given symmetric matrix.

  Parameters:
    `matrix` - The symmetric matrix to decompose.
    `splitTolerance` - tolerance on the off-diagonal elements relative to the geometric mean to split the tridiagonal matrix (a suggested value is `MathUtils.SAFE_MIN`)
  Throws:
    `InvalidMatrixException` - (wrapping a `ConvergenceException` if algorithm fails to converge

### EigenDecompositionImpl

    public EigenDecompositionImpl(double[] main,
                                  double[] secondary,
                                  double splitTolerance)
                           throws InvalidMatrixException

  Calculates the eigen decomposition of the given tridiagonal symmetric matrix.

  Parameters:
    `main` - the main diagonal of the matrix (will be copied)
    `secondary` - the secondary diagonal of the matrix (will be copied)
    `splitTolerance` - tolerance on the off-diagonal elements relative to the geometric mean to split the tridiagonal matrix (a suggested value is `MathUtils.SAFE_MIN`)
  Throws:
    `InvalidMatrixException` - (wrapping a `ConvergenceException` if algorithm fails to converge

Method Detail

### isSymmetric

    private boolean isSymmetric(RealMatrix matrix)

  Check if a matrix is symmetric.

  Parameters:
    `matrix` - matrix to check
  Returns:
    true if matrix is symmetric

### decompose

    private void decompose()

  Decompose a tridiagonal symmetric matrix.

  Throws:
    `InvalidMatrixException` - (wrapping a `ConvergenceException` if algorithm fails to converge

### getV

    public RealMatrix getV()
                    throws InvalidMatrixException

  Returns the matrix V of the decomposition.

  V is an orthogonal matrix, i.e. its transpose is also its inverse.

  The columns of V are the eigenvectors of the original matrix.

  Specified by:
    `getV` in interface `EigenDecomposition`
  Returns:
    the V matrix
  Throws:
    `InvalidMatrixException`

### getD

    public RealMatrix getD()
                    throws InvalidMatrixException

  Returns the block diagonal matrix D of the decomposition.

  D is a block diagonal matrix.

  Real eigenvalues are on the diagonal while complex values are on 2x2 blocks { {real +imaginary}, {-imaginary, real} }.

  Specified by:
    `getD` in interface `EigenDecomposition`
  Returns:
    the D matrix
  Throws:
    `InvalidMatrixException`
  See Also:
    `EigenDecomposition.getRealEigenvalues()`, `EigenDecomposition.getImagEigenvalues()`

### getVT

    public RealMatrix getVT()
                     throws InvalidMatrixException

  Returns the transpose of the matrix V of the decomposition.

  V is an orthogonal matrix, i.e. its transpose is also its inverse.

  The columns of V are the eigenvectors of the original matrix.

  Specified by:
    `getVT` in interface `EigenDecomposition`
  Returns:
    the transpose of the V matrix
  Throws:
    `InvalidMatrixException`

### getRealEigenvalues

    public double[] getRealEigenvalues()
                                throws InvalidMatrixException

  Returns a copy of the real parts of the eigenvalues of the original matrix.

  Specified by:
    `getRealEigenvalues` in interface `EigenDecomposition`
  Returns:
    a copy of the real parts of the eigenvalues of the original matrix
  Throws:
    `InvalidMatrixException`
  See Also:
    `EigenDecomposition.getD()`, `EigenDecomposition.getRealEigenvalue(int)`, `EigenDecomposition.getImagEigenvalues()`

### getRealEigenvalue

    public double getRealEigenvalue(int i)
                             throws InvalidMatrixException,
                                    java.lang.ArrayIndexOutOfBoundsException

  Returns the real part of the ith eigenvalue of the original matrix.

  Specified by:
    `getRealEigenvalue` in interface `EigenDecomposition`
  Parameters:
    `i` - index of the eigenvalue (counting from 0)
  Returns:
    real part of the ith eigenvalue of the original matrix
  Throws:
    `InvalidMatrixException`
    `java.lang.ArrayIndexOutOfBoundsException`
  See Also:
    `EigenDecomposition.getD()`, `EigenDecomposition.getRealEigenvalues()`, `EigenDecomposition.getImagEigenvalue(int)`

### getImagEigenvalues

    public double[] getImagEigenvalues()
                                throws InvalidMatrixException

  Returns a copy of the imaginary parts of the eigenvalues of the original matrix.

  Specified by:
    `getImagEigenvalues` in interface `EigenDecomposition`
  Returns:
    a copy of the imaginary parts of the eigenvalues of the original matrix
  Throws:
    `InvalidMatrixException`
  See Also:
    `EigenDecomposition.getD()`, `EigenDecomposition.getImagEigenvalue(int)`, `EigenDecomposition.getRealEigenvalues()`

### getImagEigenvalue

    public double getImagEigenvalue(int i)
                             throws InvalidMatrixException,
                                    java.lang.ArrayIndexOutOfBoundsException

  Returns the imaginary part of the ith eigenvalue of the original matrix.

  Specified by:
    `getImagEigenvalue` in interface `EigenDecomposition`
  Parameters:
    `i` - index of the eigenvalue (counting from 0)
  Returns:
    imaginary part of the ith eigenvalue of the original matrix
  Throws:
    `InvalidMatrixException`
    `java.lang.ArrayIndexOutOfBoundsException`
  See Also:
    `EigenDecomposition.getD()`, `EigenDecomposition.getImagEigenvalues()`, `EigenDecomposition.getRealEigenvalue(int)`

### getEigenvector

    public RealVector getEigenvector(int i)
                              throws InvalidMatrixException,
                                     java.lang.ArrayIndexOutOfBoundsException

  Returns a copy of the ith eigenvector of the original matrix.

  Specified by:
    `getEigenvector` in interface `EigenDecomposition`
  Parameters:
    `i` - index of the eigenvector (counting from 0)
  Returns:
    copy of the ith eigenvector of the original matrix
  Throws:
    `InvalidMatrixException`
    `java.lang.ArrayIndexOutOfBoundsException`
  See Also:
    `EigenDecomposition.getD()`

### getDeterminant

    public double getDeterminant()

  Return the determinant of the matrix

  Specified by:
    `getDeterminant` in interface `EigenDecomposition`
  Returns:
    determinant of the matrix

### getSolver

    public DecompositionSolver getSolver()

  Get a solver for finding the A × X = B solution in exact linear sense.

  Specified by:
    `getSolver` in interface `EigenDecomposition`
  Returns:
    a solver

### transformToTridiagonal

    private void transformToTridiagonal(RealMatrix matrix)

  Transform matrix to tridiagonal.

  Parameters:
    `matrix` - matrix to transform

### computeGershgorinCircles

    private void computeGershgorinCircles()

  Compute the Gershgorin circles for all rows.

### findEigenvalues

    private void findEigenvalues()
                          throws InvalidMatrixException

  Find the realEigenvalues.

  Throws:
    `InvalidMatrixException` - if a block cannot be diagonalized

### computeSplits

    private java.util.List<java.lang.Integer> computeSplits()

  Compute splitting points.

  Returns:
    list of indices after matrix can be split

### process1RowBlock

    private void process1RowBlock(int index)

  Find eigenvalue in a block with 1 row.

  In low dimensions, we simply solve the characteristic polynomial.

  Parameters:
    `index` - index of the first row of the block

### process2RowsBlock

    private void process2RowsBlock(int index)
                            throws InvalidMatrixException

  Find realEigenvalues in a block with 2 rows.

  In low dimensions, we simply solve the characteristic polynomial.

  Parameters:
    `index` - index of the first row of the block
  Throws:
    `InvalidMatrixException` - if characteristic polynomial cannot be solved

### process3RowsBlock

    private void process3RowsBlock(int index)
                            throws InvalidMatrixException

  Find realEigenvalues in a block with 3 rows.

  In low dimensions, we simply solve the characteristic polynomial.

  Parameters:
    `index` - index of the first row of the block
  Throws:
    `InvalidMatrixException` - if diagonal elements are not positive

### processGeneralBlock

    private void processGeneralBlock(int n)
                              throws InvalidMatrixException

  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.

  Parameters:
    `n` - number of rows of the block
  Throws:
    `InvalidMatrixException` - if block cannot be diagonalized after 30 * n iterations

### initialSplits

    private void initialSplits(int n)

  Perform two iterations with Li's tests for initial splits.

  Parameters:
    `n` - number of rows of the matrix to process

### goodStep

    private int goodStep(int start,
                         int end)

  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.

  Parameters:
    `start` - start index
    `end` - end index
  Returns:
    new end (maybe deflated)

### flipIfWarranted

    private boolean flipIfWarranted(int n,
                                    int step)

  Flip qd array if warranted.

  Parameters:
    `n` - number of rows in the block
    `step` - within the array (1 for flipping all elements, 2 for flipping only every other element)
  Returns:
    true if qd array was flipped

### eigenvaluesRange

    private double[] eigenvaluesRange(int index,
                                      int n)

  Compute an interval containing all realEigenvalues of a block.

  Parameters:
    `index` - index of the first row of the block
    `n` - number of rows of the block
  Returns:
    an interval containing the realEigenvalues

### countEigenValues

    private int countEigenValues(double t,
                                 int index,
                                 int n)

  Count the number of realEigenvalues below a point.

  Parameters:
    `t` - value below which we must count the number of realEigenvalues
    `index` - index of the first row of the block
    `n` - number of rows of the block
  Returns:
    number of realEigenvalues smaller than t

### ldlTDecomposition

    private void ldlTDecomposition(double lambda,
                                   int index,
                                   int n)

  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.

  Parameters:
    `lambda` - shift to add to the matrix before decomposing it to ensure it is positive definite
    `index` - index of the first row of the block
    `n` - number of rows of the block

### dqds

    private void dqds(int start,
                      int end)

  Perform a dqds step, using current shift increment.

  This implementation is a translation of the LAPACK routine DLASQ5.

  Parameters:
    `start` - start index
    `end` - end index

### dqd

    private void dqd(int start,
                     int end)

  Perform a dqd step.

  This implementation is a translation of the LAPACK routine DLASQ6.

  Parameters:
    `start` - start index
    `end` - end index

### computeShiftIncrement

    private void computeShiftIncrement(int start,
                                       int end,
                                       int deflated)

  Compute the shift increment as an estimate of the smallest eigenvalue.

  This implementation is a translation of the LAPACK routine DLAZQ4.

  Parameters:
    `start` - start index
    `end` - end index
    `deflated` - number of realEigenvalues just deflated

### updateSigma

    private void updateSigma(double tau)

  Update sigma.

  Parameters:
    `tau` - shift to apply to sigma

### findEigenVectors

    private void findEigenVectors()

  Find eigenvectors.

### findEigenvector

    private ArrayRealVector findEigenvector(double eigenvalue,
                                            double[] d,
                                            double[] l)

  Find an eigenvector corresponding to an eigenvalue, using bidiagonals.

  This method corresponds to algorithm X from Dhillon's thesis.

  Parameters:
    `eigenvalue` - eigenvalue for which eigenvector is desired
    `d` - diagonal elements of the initial non-shifted D matrix
    `l` - off-diagonal elements of the initial non-shifted L matrix
  Returns:
    an eigenvector

### stationaryQuotientDifferenceWithShift

    private void stationaryQuotientDifferenceWithShift(double[] d,
                                                       double[] l,
                                                       double lambda)

  Decompose matrix LDLT - λ I as L+D+L+T.

  This method corresponds to algorithm 4.4.3 (dstqds) from Dhillon's thesis.

  Parameters:
    `d` - diagonal elements of D,
    `l` - off-diagonal elements of L
    `lambda` - shift to apply

### progressiveQuotientDifferenceWithShift

    private void progressiveQuotientDifferenceWithShift(double[] d,
                                                        double[] l,
                                                        double lambda)

  Decompose matrix LDLT - λ I as U-D-U-T.

  This method corresponds to algorithm 4.4.5 (dqds) from Dhillon's thesis.

  Parameters:
    `d` - diagonal elements of D
    `l` - off-diagonal elements of L
    `lambda` - shift to apply

Copyright (c) 2003-2010 Apache Software Foundation