org.apache.commons.math.linear
Class EigenDecompositionImpl

java.lang.Object
  extended by org.apache.commons.math.linear.EigenDecompositionImpl
All Implemented Interfaces:
EigenDecomposition

public class EigenDecompositionImpl
extends 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  List<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,
                                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
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,
                                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
ArrayIndexOutOfBoundsException
See Also:
EigenDecomposition.getD(), EigenDecomposition.getImagEigenvalues(), EigenDecomposition.getRealEigenvalue(int)

getEigenvector

public RealVector getEigenvector(int i)
                          throws InvalidMatrixException,
                                 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
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 List<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