src/main/java/org/apache/commons/math3/linear/EigenDecomposition.java - commons-math - Git at Google

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

 package org.apache.commons.math3.linear;

 import org.apache.commons.math3.exception.MaxCountExceededException;
 import org.apache.commons.math3.exception.DimensionMismatchException;
 import org.apache.commons.math3.exception.util.LocalizedFormats;
 import org.apache.commons.math3.util.Precision;
 import org.apache.commons.math3.util.FastMath;

 /**
  * Calculates the eigen decomposition of a real <strong>symmetric</strong>
  * matrix.
  * <p>The eigen decomposition of matrix A is a set of two matrices:
  * V and D such that A = V &times; D &times; V<sup>T</sup>.
  * A, V and D are all m &times; m matrices.</p>
  * <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code>
  * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
  * library, with the following changes:</p>
  * <ul>
  *   <li>a {@link #getVT() getVt} method has been added,</li>
  *   <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int)
  *   getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li>
  *   <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single
  *   eigenvector has been added,</li>
  *   <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
  *   <li>a {@link #getSolver() getSolver} method has been added.</li>
  * </ul>
  * <p>
  * As of 2.0, this class supports only <strong>symmetric</strong> matrices, and
  * hence computes only real realEigenvalues. This implies the D matrix returned
  * by {@link #getD()} is always diagonal and the imaginary values returned
  * {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always
  * null.
  * </p>
  * <p>
  * When called with a {@link RealMatrix} argument, this implementation only uses
  * the upper part of the matrix, the part below the diagonal is not accessed at
  * all.
  * </p>
  * <p>
  * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
  * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
  * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
  * New-York
  * </p>
  * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
  * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
  * @version $Id$
  * @since 2.0 (changed to concrete class in 3.0)
  */
 public class EigenDecomposition{
     /** Maximum number of iterations accepted in the implicit QL transformation */
     private byte maxIter = 30;
     /** Main diagonal of the tridiagonal matrix. */
     private double[] main;
     /** Secondary diagonal of the tridiagonal matrix. */
     private double[] secondary;
     /**
      * Transformer to tridiagonal (may be null if matrix is already
      * tridiagonal).
      */
     private TriDiagonalTransformer transformer;
     /** Real part of the realEigenvalues. */
     private double[] realEigenvalues;
     /** Imaginary part of the realEigenvalues. */
     private double[] imagEigenvalues;
     /** Eigenvectors. */
     private ArrayRealVector[] eigenvectors;
     /** Cached value of V. */
     private RealMatrix cachedV;
     /** Cached value of D. */
     private RealMatrix cachedD;
     /** Cached value of Vt. */
     private RealMatrix cachedVt;

     /**
      * Calculates the eigen decomposition of the given symmetric matrix.
      *
      * @param matrix Matrix to decompose. It <em>must</em> be symmetric.
      * @param splitTolerance Dummy parameter (present for backward
      * compatibility only).
      * @throws NonSymmetricMatrixException if the matrix is not symmetric.
      * @throws MaxCountExceededException if the algorithm fails to converge.
      */
     public EigenDecomposition(final RealMatrix matrix,
                                   final double splitTolerance)  {
         if (isSymmetric(matrix, true)) {
             transformToTridiagonal(matrix);
             findEigenVectors(transformer.getQ().getData());
         }
     }

     /**
      * Calculates the eigen decomposition of the symmetric tridiagonal
      * matrix.  The Householder matrix is assumed to be the identity matrix.
      *
      * @param main Main diagonal of the symmetric tridiagonal form.
      * @param secondary Secondary of the tridiagonal form.
      * @param splitTolerance Dummy parameter (present for backward
      * compatibility only).
      * @throws MaxCountExceededException if the algorithm fails to converge.
      */
     public EigenDecomposition(final double[] main,final double[] secondary,
                                   final double splitTolerance) {
         this.main      = main.clone();
         this.secondary = secondary.clone();
         transformer    = null;
         final int size=main.length;
         double[][] z = new double[size][size];
         for (int i=0;i<size;i++) {
             z[i][i]=1.0;
         }
         findEigenVectors(z);
     }

     /**
      * Check if a matrix is symmetric.
      *
      * @param matrix Matrix to check.
      * @param raiseException If {@code true}, the method will throw an
      * exception if {@code matrix} is not symmetric.
      * @return {@code true} if {@code matrix} is symmetric.
      * @throws NonSymmetricMatrixException if the matrix is not symmetric and
      * {@code raiseException} is {@code true}.
      */
     private boolean isSymmetric(final RealMatrix matrix,
                                 boolean raiseException) {
         final int rows = matrix.getRowDimension();
         final int columns = matrix.getColumnDimension();
         final double eps = 10 * rows * columns * Precision.EPSILON;
         for (int i = 0; i < rows; ++i) {
             for (int j = i + 1; j < columns; ++j) {
                 final double mij = matrix.getEntry(i, j);
                 final double mji = matrix.getEntry(j, i);
                 if (FastMath.abs(mij - mji) >
                     (FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
                     if (raiseException) {
                         throw new NonSymmetricMatrixException(i, j, eps);
                     }
                     return false;
                 }
             }
         }
         return true;
     }

     /**
      * Gets 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.
      * No assumption is made about the orientation of the system axes formed
      * by the columns of V (e.g. in a 3-dimension space, V can form a left-
      * or right-handed system).
      *
      * @return the V matrix.
      */
     public RealMatrix getV() {

         if (cachedV == null) {
             final int m = eigenvectors.length;
             cachedV = MatrixUtils.createRealMatrix(m, m);
             for (int k = 0; k < m; ++k) {
                 cachedV.setColumnVector(k, eigenvectors[k]);
             }
         }
         // return the cached matrix
         return cachedV;

     }

     /**
      * Gets 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} }.
      *
      * @return the D matrix.
      *
      * @see #getRealEigenvalues()
      * @see #getImagEigenvalues()
      */
     public RealMatrix getD() {
         if (cachedD == null) {
             // cache the matrix for subsequent calls
             cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
         }
         return cachedD;
     }

     /**
      * Gets 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.
      * No assumption is made about the orientation of the system axes formed
      * by the columns of V (e.g. in a 3-dimension space, V can form a left-
      * or right-handed system).
      *
      * @return the transpose of the V matrix.
      */
     public RealMatrix getVT() {

         if (cachedVt == null) {
             final int m = eigenvectors.length;
             cachedVt = MatrixUtils.createRealMatrix(m, m);
             for (int k = 0; k < m; ++k) {
                 cachedVt.setRowVector(k, eigenvectors[k]);
             }

         }

         // return the cached matrix
         return cachedVt;
     }

     /**
      * Gets a copy of the real parts of the eigenvalues of the original matrix.
      *
      * @return a copy of the real parts of the eigenvalues of the original matrix.
      *
      * @see #getD()
      * @see #getRealEigenvalue(int)
      * @see #getImagEigenvalues()
      */
     public double[] getRealEigenvalues() {
         return realEigenvalues.clone();
     }

     /**
      * Returns the real part of the i<sup>th</sup> eigenvalue of the original
      * matrix.
      *
      * @param i index of the eigenvalue (counting from 0)
      * @return real part of the i<sup>th</sup> eigenvalue of the original
      * matrix.
      *
      * @see #getD()
      * @see #getRealEigenvalues()
      * @see #getImagEigenvalue(int)
      */
     public double getRealEigenvalue(final int i) {
         return realEigenvalues[i];
     }

     /**
      * Gets a copy of the imaginary parts of the eigenvalues of the original
      * matrix.
      *
      * @return a copy of the imaginary parts of the eigenvalues of the original
      * matrix.
      *
      * @see #getD()
      * @see #getImagEigenvalue(int)
      * @see #getRealEigenvalues()
      */
     public double[] getImagEigenvalues() {
         return imagEigenvalues.clone();
     }

     /**
      * Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original
      * matrix.
      *
      * @param i Index of the eigenvalue (counting from 0).
      * @return the imaginary part of the i<sup>th</sup> eigenvalue of the original
      * matrix.
      *
      * @see #getD()
      * @see #getImagEigenvalues()
      * @see #getRealEigenvalue(int)
      */
     public double getImagEigenvalue(final int i) {
         return imagEigenvalues[i];
     }

     /**
      * Gets a copy of the i<sup>th</sup> eigenvector of the original matrix.
      *
      * @param i Index of the eigenvector (counting from 0).
      * @return a copy of the i<sup>th</sup> eigenvector of the original matrix.
      * @see #getD()
      */
     public RealVector getEigenvector(final int i) {
         return eigenvectors[i].copy();
     }

     /**
      * Computes the determinant of the matrix.
      *
      * @return the determinant of the matrix.
      */
     public double getDeterminant() {
         double determinant = 1;
         for (double lambda : realEigenvalues) {
             determinant *= lambda;
         }
         return determinant;
     }

     /**
      * Gets a solver for finding the A &times; X = B solution in exact
      * linear sense.
      *
      * @return a solver.
      */
     public DecompositionSolver getSolver() {
         return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
     }

     /** Specialized solver. */
     private static class Solver implements DecompositionSolver {
         /** Real part of the realEigenvalues. */
         private double[] realEigenvalues;
         /** Imaginary part of the realEigenvalues. */
         private double[] imagEigenvalues;
         /** Eigenvectors. */
         private final ArrayRealVector[] eigenvectors;

         /**
          * Builds a solver from decomposed matrix.
          *
          * @param realEigenvalues Real parts of the eigenvalues.
          * @param imagEigenvalues Imaginary parts of the eigenvalues.
          * @param eigenvectors Eigenvectors.
          */
         private Solver(final double[] realEigenvalues,
                 final double[] imagEigenvalues,
                 final ArrayRealVector[] eigenvectors) {
             this.realEigenvalues = realEigenvalues;
             this.imagEigenvalues = imagEigenvalues;
             this.eigenvectors = eigenvectors;
         }

         /**
          * Solves the linear equation A &times; X = B for symmetric matrices A.
          * <p>
          * This method only finds exact linear solutions, i.e. solutions for
          * which ||A &times; X - B|| is exactly 0.
          * </p>
          *
          * @param b Right-hand side of the equation A &times; X = B.
          * @return a Vector X that minimizes the two norm of A &times; X - B.
          *
          * @throws DimensionMismatchException if the matrices dimensions do not match.
          * @throws SingularMatrixException if the decomposed matrix is singular.
          */
         public RealVector solve(final RealVector b) {
             if (!isNonSingular()) {
                 throw new SingularMatrixException();
             }

             final int m = realEigenvalues.length;
             if (b.getDimension() != m) {
                 throw new DimensionMismatchException(b.getDimension(), m);
             }

             final double[] bp = new double[m];
             for (int i = 0; i < m; ++i) {
                 final ArrayRealVector v = eigenvectors[i];
                 final double[] vData = v.getDataRef();
                 final double s = v.dotProduct(b) / realEigenvalues[i];
                 for (int j = 0; j < m; ++j) {
                     bp[j] += s * vData[j];
                 }
             }

             return new ArrayRealVector(bp, false);
         }

         /** {@inheritDoc} */
         public RealMatrix solve(RealMatrix b) {

             if (!isNonSingular()) {
                 throw new SingularMatrixException();
             }

             final int m = realEigenvalues.length;
             if (b.getRowDimension() != m) {
                 throw new DimensionMismatchException(b.getRowDimension(), m);
             }

             final int nColB = b.getColumnDimension();
             final double[][] bp = new double[m][nColB];
             final double[] tmpCol = new double[m];
             for (int k = 0; k < nColB; ++k) {
                 for (int i = 0; i < m; ++i) {
                     tmpCol[i] = b.getEntry(i, k);
                     bp[i][k]  = 0;
                 }
                 for (int i = 0; i < m; ++i) {
                     final ArrayRealVector v = eigenvectors[i];
                     final double[] vData = v.getDataRef();
                     double s = 0;
                     for (int j = 0; j < m; ++j) {
                         s += v.getEntry(j) * tmpCol[j];
                     }
                     s /= realEigenvalues[i];
                     for (int j = 0; j < m; ++j) {
                         bp[j][k] += s * vData[j];
                     }
                 }
             }

             return new Array2DRowRealMatrix(bp, false);

         }

         /**
          * Checks whether the decomposed matrix is non-singular.
          *
          * @return true if the decomposed matrix is non-singular.
          */
         public boolean isNonSingular() {
             for (int i = 0; i < realEigenvalues.length; ++i) {
                 if (realEigenvalues[i] == 0 &&
                     imagEigenvalues[i] == 0) {
                     return false;
                 }
             }
             return true;
         }

         /**
          * Get the inverse of the decomposed matrix.
          *
          * @return the inverse matrix.
          * @throws SingularMatrixException if the decomposed matrix is singular.
          */
         public RealMatrix getInverse() {
             if (!isNonSingular()) {
                 throw new SingularMatrixException();
             }

             final int m = realEigenvalues.length;
             final double[][] invData = new double[m][m];

             for (int i = 0; i < m; ++i) {
                 final double[] invI = invData[i];
                 for (int j = 0; j < m; ++j) {
                     double invIJ = 0;
                     for (int k = 0; k < m; ++k) {
                         final double[] vK = eigenvectors[k].getDataRef();
                         invIJ += vK[i] * vK[j] / realEigenvalues[k];
                     }
                     invI[j] = invIJ;
                 }
             }
             return MatrixUtils.createRealMatrix(invData);
         }
     }

     /**
      * Transforms the matrix to tridiagonal form.
      *
      * @param matrix Matrix to transform.
      */
     private void transformToTridiagonal(final RealMatrix matrix) {
         // transform the matrix to tridiagonal
         transformer = new TriDiagonalTransformer(matrix);
         main = transformer.getMainDiagonalRef();
         secondary = transformer.getSecondaryDiagonalRef();
     }

     /**
      * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
      *
      * @param householderMatrix Householder matrix of the transformation
      * to tridiagonal form.
      */
     private void findEigenVectors(double[][] householderMatrix) {
         final double[][]z = householderMatrix.clone();
         final int n = main.length;
         realEigenvalues = new double[n];
         imagEigenvalues = new double[n];
         final double[] e = new double[n];
         for (int i = 0; i < n - 1; i++) {
             realEigenvalues[i] = main[i];
             e[i] = secondary[i];
         }
         realEigenvalues[n - 1] = main[n - 1];
         e[n - 1] = 0;

         // Determine the largest main and secondary value in absolute term.
         double maxAbsoluteValue = 0;
         for (int i = 0; i < n; i++) {
             if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
                 maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
             }
             if (FastMath.abs(e[i]) > maxAbsoluteValue) {
                 maxAbsoluteValue = FastMath.abs(e[i]);
             }
         }
         // Make null any main and secondary value too small to be significant
         if (maxAbsoluteValue != 0) {
             for (int i=0; i < n; i++) {
                 if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
                     realEigenvalues[i] = 0;
                 }
                 if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
                     e[i]=0;
                 }
             }
         }

         for (int j = 0; j < n; j++) {
             int its = 0;
             int m;
             do {
                 for (m = j; m < n - 1; m++) {
                     double delta = FastMath.abs(realEigenvalues[m]) +
                         FastMath.abs(realEigenvalues[m + 1]);
                     if (FastMath.abs(e[m]) + delta == delta) {
                         break;
                     }
                 }
                 if (m != j) {
                     if (its == maxIter) {
                         throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
                                                             maxIter);
                     }
                     its++;
                     double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
                     double t = FastMath.sqrt(1 + q * q);
                     if (q < 0.0) {
                         q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
                     } else {
                         q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
                     }
                     double u = 0.0;
                     double s = 1.0;
                     double c = 1.0;
                     int i;
                     for (i = m - 1; i >= j; i--) {
                         double p = s * e[i];
                         double h = c * e[i];
                         if (FastMath.abs(p) >= FastMath.abs(q)) {
                             c = q / p;
                             t = FastMath.sqrt(c * c + 1.0);
                             e[i + 1] = p * t;
                             s = 1.0 / t;
                             c = c * s;
                         } else {
                             s = p / q;
                             t = FastMath.sqrt(s * s + 1.0);
                             e[i + 1] = q * t;
                             c = 1.0 / t;
                             s = s * c;
                         }
                         if (e[i + 1] == 0.0) {
                             realEigenvalues[i + 1] -= u;
                             e[m] = 0.0;
                             break;
                         }
                         q = realEigenvalues[i + 1] - u;
                         t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
                         u = s * t;
                         realEigenvalues[i + 1] = q + u;
                         q = c * t - h;
                         for (int ia = 0; ia < n; ia++) {
                             p = z[ia][i + 1];
                             z[ia][i + 1] = s * z[ia][i] + c * p;
                             z[ia][i] = c * z[ia][i] - s * p;
                         }
                     }
                     if (t == 0.0 && i >= j) {
                         continue;
                     }
                     realEigenvalues[j] -= u;
                     e[j] = q;
                     e[m] = 0.0;
                 }
             } while (m != j);
         }

         //Sort the eigen values (and vectors) in increase order
         for (int i = 0; i < n; i++) {
             int k = i;
             double p = realEigenvalues[i];
             for (int j = i + 1; j < n; j++) {
                 if (realEigenvalues[j] > p) {
                     k = j;
                     p = realEigenvalues[j];
                 }
             }
             if (k != i) {
                 realEigenvalues[k] = realEigenvalues[i];
                 realEigenvalues[i] = p;
                 for (int j = 0; j < n; j++) {
                     p = z[j][i];
                     z[j][i] = z[j][k];
                     z[j][k] = p;
                 }
             }
         }

         // Determine the largest eigen value in absolute term.
         maxAbsoluteValue = 0;
         for (int i = 0; i < n; i++) {
             if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
                 maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
             }
         }
         // Make null any eigen value too small to be significant
         if (maxAbsoluteValue!=0.0) {
             for (int i=0; i < n; i++) {
                 if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
                     realEigenvalues[i] = 0;
                 }
             }
         }
         eigenvectors = new ArrayRealVector[n];
         final double[] tmp = new double[n];
         for (int i = 0; i < n; i++) {
             for (int j = 0; j < n; j++) {
                 tmp[j] = z[j][i];
             }
             eigenvectors[i] = new ArrayRealVector(tmp);
         }
     }
 }
	/*
	* Licensed to the Apache Software Foundation (ASF) under one or more
	* contributor license agreements. See the NOTICE file distributed with
	* this work for additional information regarding copyright ownership.
	* The ASF licenses this file to You under the Apache License, Version 2.0
	* (the "License"); you may not use this file except in compliance with
	* the License. You may obtain a copy of the License at
	*
	* http://www.apache.org/licenses/LICENSE-2.0
	*
	* Unless required by applicable law or agreed to in writing, software
	* distributed under the License is distributed on an "AS IS" BASIS,
	* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	* See the License for the specific language governing permissions and
	* limitations under the License.
	*/

	package org.apache.commons.math3.linear;

	import org.apache.commons.math3.exception.MaxCountExceededException;
	import org.apache.commons.math3.exception.DimensionMismatchException;
	import org.apache.commons.math3.exception.util.LocalizedFormats;
	import org.apache.commons.math3.util.Precision;
	import org.apache.commons.math3.util.FastMath;

	/**
	* Calculates the eigen decomposition of a real <strong>symmetric</strong>
	* matrix.
	* <p>The eigen decomposition of matrix A is a set of two matrices:
	* V and D such that A = V × D × V<sup>T</sup>.
	* A, V and D are all m × m matrices.</p>
	* <p>This class is similar in spirit to the <code>EigenvalueDecomposition</code>
	* class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
	* library, with the following changes:</p>
	* <ul>
	* <li>a {@link #getVT() getVt} method has been added,</li>
	* <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and {@link #getImagEigenvalue(int)
	* getImagEigenvalue} methods to pick up a single eigenvalue have been added,</li>
	* <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a single
	* eigenvector has been added,</li>
	* <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
	* <li>a {@link #getSolver() getSolver} method has been added.</li>
	* </ul>
	* <p>
	* As of 2.0, this class supports only <strong>symmetric</strong> matrices, and
	* hence computes only real realEigenvalues. This implies the D matrix returned
	* by {@link #getD()} is always diagonal and the imaginary values returned
	* {@link #getImagEigenvalue(int)} and {@link #getImagEigenvalues()} are always
	* null.
	* </p>
	* <p>
	* When called with a {@link RealMatrix} argument, this implementation only uses
	* the upper part of the matrix, the part below the diagonal is not accessed at
	* all.
	* </p>
	* <p>
	* This implementation is based on the paper by A. Drubrulle, R.S. Martin and
	* J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
	* Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
	* New-York
	* </p>
	* @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
	* @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
	* @version $Id$
	* @since 2.0 (changed to concrete class in 3.0)
	*/
	public class EigenDecomposition{
	/** Maximum number of iterations accepted in the implicit QL transformation */
	private byte maxIter = 30;
	/** Main diagonal of the tridiagonal matrix. */
	private double[] main;
	/** Secondary diagonal of the tridiagonal matrix. */
	private double[] secondary;
	/**
	* Transformer to tridiagonal (may be null if matrix is already
	* tridiagonal).
	*/
	private TriDiagonalTransformer transformer;
	/** Real part of the realEigenvalues. */
	private double[] realEigenvalues;
	/** Imaginary part of the realEigenvalues. */
	private double[] imagEigenvalues;
	/** Eigenvectors. */
	private ArrayRealVector[] eigenvectors;
	/** Cached value of V. */
	private RealMatrix cachedV;
	/** Cached value of D. */
	private RealMatrix cachedD;
	/** Cached value of Vt. */
	private RealMatrix cachedVt;

	/**
	* Calculates the eigen decomposition of the given symmetric matrix.
	*
	* @param matrix Matrix to decompose. It <em>must</em> be symmetric.
	* @param splitTolerance Dummy parameter (present for backward
	* compatibility only).
	* @throws NonSymmetricMatrixException if the matrix is not symmetric.
	* @throws MaxCountExceededException if the algorithm fails to converge.
	*/
	public EigenDecomposition(final RealMatrix matrix,
	final double splitTolerance) {
	if (isSymmetric(matrix, true)) {
	transformToTridiagonal(matrix);
	findEigenVectors(transformer.getQ().getData());
	}
	}

	/**
	* Calculates the eigen decomposition of the symmetric tridiagonal
	* matrix. The Householder matrix is assumed to be the identity matrix.
	*
	* @param main Main diagonal of the symmetric tridiagonal form.
	* @param secondary Secondary of the tridiagonal form.
	* @param splitTolerance Dummy parameter (present for backward
	* compatibility only).
	* @throws MaxCountExceededException if the algorithm fails to converge.
	*/
	public EigenDecomposition(final double[] main,final double[] secondary,
	final double splitTolerance) {
	this.main = main.clone();
	this.secondary = secondary.clone();
	transformer = null;
	final int size=main.length;
	double[][] z = new double[size][size];
	for (int i=0;i<size;i++) {
	z[i][i]=1.0;
	}
	findEigenVectors(z);
	}

	/**
	* Check if a matrix is symmetric.
	*
	* @param matrix Matrix to check.
	* @param raiseException If {@code true}, the method will throw an
	* exception if {@code matrix} is not symmetric.
	* @return {@code true} if {@code matrix} is symmetric.
	* @throws NonSymmetricMatrixException if the matrix is not symmetric and
	* {@code raiseException} is {@code true}.
	*/
	private boolean isSymmetric(final RealMatrix matrix,
	boolean raiseException) {
	final int rows = matrix.getRowDimension();
	final int columns = matrix.getColumnDimension();
	final double eps = 10 * rows * columns * Precision.EPSILON;
	for (int i = 0; i < rows; ++i) {
	for (int j = i + 1; j < columns; ++j) {
	final double mij = matrix.getEntry(i, j);
	final double mji = matrix.getEntry(j, i);
	if (FastMath.abs(mij - mji) >
	(FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
	if (raiseException) {
	throw new NonSymmetricMatrixException(i, j, eps);
	}
	return false;
	}
	}
	}
	return true;
	}

	/**
	* Gets 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.
	* No assumption is made about the orientation of the system axes formed
	* by the columns of V (e.g. in a 3-dimension space, V can form a left-
	* or right-handed system).
	*
	* @return the V matrix.
	*/
	public RealMatrix getV() {

	if (cachedV == null) {
	final int m = eigenvectors.length;
	cachedV = MatrixUtils.createRealMatrix(m, m);
	for (int k = 0; k < m; ++k) {
	cachedV.setColumnVector(k, eigenvectors[k]);
	}
	}
	// return the cached matrix
	return cachedV;

	}

	/**
	* Gets 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} }.
	*
	* @return the D matrix.
	*
	* @see #getRealEigenvalues()
	* @see #getImagEigenvalues()
	*/
	public RealMatrix getD() {
	if (cachedD == null) {
	// cache the matrix for subsequent calls
	cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
	}
	return cachedD;
	}

	/**
	* Gets 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.
	* No assumption is made about the orientation of the system axes formed
	* by the columns of V (e.g. in a 3-dimension space, V can form a left-
	* or right-handed system).
	*
	* @return the transpose of the V matrix.
	*/
	public RealMatrix getVT() {

	if (cachedVt == null) {
	final int m = eigenvectors.length;
	cachedVt = MatrixUtils.createRealMatrix(m, m);
	for (int k = 0; k < m; ++k) {
	cachedVt.setRowVector(k, eigenvectors[k]);
	}

	}

	// return the cached matrix
	return cachedVt;
	}

	/**
	* Gets a copy of the real parts of the eigenvalues of the original matrix.
	*
	* @return a copy of the real parts of the eigenvalues of the original matrix.
	*
	* @see #getD()
	* @see #getRealEigenvalue(int)
	* @see #getImagEigenvalues()
	*/
	public double[] getRealEigenvalues() {
	return realEigenvalues.clone();
	}

	/**
	* Returns the real part of the i<sup>th</sup> eigenvalue of the original
	* matrix.
	*
	* @param i index of the eigenvalue (counting from 0)
	* @return real part of the i<sup>th</sup> eigenvalue of the original
	* matrix.
	*
	* @see #getD()
	* @see #getRealEigenvalues()
	* @see #getImagEigenvalue(int)
	*/
	public double getRealEigenvalue(final int i) {
	return realEigenvalues[i];
	}

	/**
	* Gets a copy of the imaginary parts of the eigenvalues of the original
	* matrix.
	*
	* @return a copy of the imaginary parts of the eigenvalues of the original
	* matrix.
	*
	* @see #getD()
	* @see #getImagEigenvalue(int)
	* @see #getRealEigenvalues()
	*/
	public double[] getImagEigenvalues() {
	return imagEigenvalues.clone();
	}

	/**
	* Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original
	* matrix.
	*
	* @param i Index of the eigenvalue (counting from 0).
	* @return the imaginary part of the i<sup>th</sup> eigenvalue of the original
	* matrix.
	*
	* @see #getD()
	* @see #getImagEigenvalues()
	* @see #getRealEigenvalue(int)
	*/
	public double getImagEigenvalue(final int i) {
	return imagEigenvalues[i];
	}

	/**
	* Gets a copy of the i<sup>th</sup> eigenvector of the original matrix.
	*
	* @param i Index of the eigenvector (counting from 0).
	* @return a copy of the i<sup>th</sup> eigenvector of the original matrix.
	* @see #getD()
	*/
	public RealVector getEigenvector(final int i) {
	return eigenvectors[i].copy();
	}

	/**
	* Computes the determinant of the matrix.
	*
	* @return the determinant of the matrix.
	*/
	public double getDeterminant() {
	double determinant = 1;
	for (double lambda : realEigenvalues) {
	determinant *= lambda;
	}
	return determinant;
	}

	/**
	* Gets a solver for finding the A × X = B solution in exact
	* linear sense.
	*
	* @return a solver.
	*/
	public DecompositionSolver getSolver() {
	return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
	}

	/** Specialized solver. */
	private static class Solver implements DecompositionSolver {
	/** Real part of the realEigenvalues. */
	private double[] realEigenvalues;
	/** Imaginary part of the realEigenvalues. */
	private double[] imagEigenvalues;
	/** Eigenvectors. */
	private final ArrayRealVector[] eigenvectors;

	/**
	* Builds a solver from decomposed matrix.
	*
	* @param realEigenvalues Real parts of the eigenvalues.
	* @param imagEigenvalues Imaginary parts of the eigenvalues.
	* @param eigenvectors Eigenvectors.
	*/
	private Solver(final double[] realEigenvalues,
	final double[] imagEigenvalues,
	final ArrayRealVector[] eigenvectors) {
	this.realEigenvalues = realEigenvalues;
	this.imagEigenvalues = imagEigenvalues;
	this.eigenvectors = eigenvectors;
	}

	/**
	* Solves the linear equation A × X = B for symmetric matrices A.
	* <p>
	* This method only finds exact linear solutions, i.e. solutions for
	* which \|\|A × X - B\|\| is exactly 0.
	* </p>
	*
	* @param b Right-hand side of the equation A × X = B.
	* @return a Vector X that minimizes the two norm of A × X - B.
	*
	* @throws DimensionMismatchException if the matrices dimensions do not match.
	* @throws SingularMatrixException if the decomposed matrix is singular.
	*/
	public RealVector solve(final RealVector b) {
	if (!isNonSingular()) {
	throw new SingularMatrixException();
	}

	final int m = realEigenvalues.length;
	if (b.getDimension() != m) {
	throw new DimensionMismatchException(b.getDimension(), m);
	}

	final double[] bp = new double[m];
	for (int i = 0; i < m; ++i) {
	final ArrayRealVector v = eigenvectors[i];
	final double[] vData = v.getDataRef();
	final double s = v.dotProduct(b) / realEigenvalues[i];
	for (int j = 0; j < m; ++j) {
	bp[j] += s * vData[j];
	}
	}

	return new ArrayRealVector(bp, false);
	}

	/** {@inheritDoc} */
	public RealMatrix solve(RealMatrix b) {

	if (!isNonSingular()) {
	throw new SingularMatrixException();
	}

	final int m = realEigenvalues.length;
	if (b.getRowDimension() != m) {
	throw new DimensionMismatchException(b.getRowDimension(), m);
	}

	final int nColB = b.getColumnDimension();
	final double[][] bp = new double[m][nColB];
	final double[] tmpCol = new double[m];
	for (int k = 0; k < nColB; ++k) {
	for (int i = 0; i < m; ++i) {
	tmpCol[i] = b.getEntry(i, k);
	bp[i][k] = 0;
	}
	for (int i = 0; i < m; ++i) {
	final ArrayRealVector v = eigenvectors[i];
	final double[] vData = v.getDataRef();
	double s = 0;
	for (int j = 0; j < m; ++j) {
	s += v.getEntry(j) * tmpCol[j];
	}
	s /= realEigenvalues[i];
	for (int j = 0; j < m; ++j) {
	bp[j][k] += s * vData[j];
	}
	}
	}

	return new Array2DRowRealMatrix(bp, false);

	}

	/**
	* Checks whether the decomposed matrix is non-singular.
	*
	* @return true if the decomposed matrix is non-singular.
	*/
	public boolean isNonSingular() {
	for (int i = 0; i < realEigenvalues.length; ++i) {
	if (realEigenvalues[i] == 0 &&
	imagEigenvalues[i] == 0) {
	return false;
	}
	}
	return true;
	}

	/**
	* Get the inverse of the decomposed matrix.
	*
	* @return the inverse matrix.
	* @throws SingularMatrixException if the decomposed matrix is singular.
	*/
	public RealMatrix getInverse() {
	if (!isNonSingular()) {
	throw new SingularMatrixException();
	}

	final int m = realEigenvalues.length;
	final double[][] invData = new double[m][m];

	for (int i = 0; i < m; ++i) {
	final double[] invI = invData[i];
	for (int j = 0; j < m; ++j) {
	double invIJ = 0;
	for (int k = 0; k < m; ++k) {
	final double[] vK = eigenvectors[k].getDataRef();
	invIJ += vK[i] * vK[j] / realEigenvalues[k];
	}
	invI[j] = invIJ;
	}
	}
	return MatrixUtils.createRealMatrix(invData);
	}
	}

	/**
	* Transforms the matrix to tridiagonal form.
	*
	* @param matrix Matrix to transform.
	*/
	private void transformToTridiagonal(final RealMatrix matrix) {
	// transform the matrix to tridiagonal
	transformer = new TriDiagonalTransformer(matrix);
	main = transformer.getMainDiagonalRef();
	secondary = transformer.getSecondaryDiagonalRef();
	}

	/**
	* Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
	*
	* @param householderMatrix Householder matrix of the transformation
	* to tridiagonal form.
	*/
	private void findEigenVectors(double[][] householderMatrix) {
	final double[][]z = householderMatrix.clone();
	final int n = main.length;
	realEigenvalues = new double[n];
	imagEigenvalues = new double[n];
	final double[] e = new double[n];
	for (int i = 0; i < n - 1; i++) {
	realEigenvalues[i] = main[i];
	e[i] = secondary[i];
	}
	realEigenvalues[n - 1] = main[n - 1];
	e[n - 1] = 0;

	// Determine the largest main and secondary value in absolute term.
	double maxAbsoluteValue = 0;
	for (int i = 0; i < n; i++) {
	if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
	maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
	}
	if (FastMath.abs(e[i]) > maxAbsoluteValue) {
	maxAbsoluteValue = FastMath.abs(e[i]);
	}
	}
	// Make null any main and secondary value too small to be significant
	if (maxAbsoluteValue != 0) {
	for (int i=0; i < n; i++) {
	if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
	realEigenvalues[i] = 0;
	}
	if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
	e[i]=0;
	}
	}
	}

	for (int j = 0; j < n; j++) {
	int its = 0;
	int m;
	do {
	for (m = j; m < n - 1; m++) {
	double delta = FastMath.abs(realEigenvalues[m]) +
	FastMath.abs(realEigenvalues[m + 1]);
	if (FastMath.abs(e[m]) + delta == delta) {
	break;
	}
	}
	if (m != j) {
	if (its == maxIter) {
	throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
	maxIter);
	}
	its++;
	double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
	double t = FastMath.sqrt(1 + q * q);
	if (q < 0.0) {
	q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
	} else {
	q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
	}
	double u = 0.0;
	double s = 1.0;
	double c = 1.0;
	int i;
	for (i = m - 1; i >= j; i--) {
	double p = s * e[i];
	double h = c * e[i];
	if (FastMath.abs(p) >= FastMath.abs(q)) {
	c = q / p;
	t = FastMath.sqrt(c * c + 1.0);
	e[i + 1] = p * t;
	s = 1.0 / t;
	c = c * s;
	} else {
	s = p / q;
	t = FastMath.sqrt(s * s + 1.0);
	e[i + 1] = q * t;
	c = 1.0 / t;
	s = s * c;
	}
	if (e[i + 1] == 0.0) {
	realEigenvalues[i + 1] -= u;
	e[m] = 0.0;
	break;
	}
	q = realEigenvalues[i + 1] - u;
	t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
	u = s * t;
	realEigenvalues[i + 1] = q + u;
	q = c * t - h;
	for (int ia = 0; ia < n; ia++) {
	p = z[ia][i + 1];
	z[ia][i + 1] = s * z[ia][i] + c * p;
	z[ia][i] = c * z[ia][i] - s * p;
	}
	}
	if (t == 0.0 && i >= j) {
	continue;
	}
	realEigenvalues[j] -= u;
	e[j] = q;
	e[m] = 0.0;
	}
	} while (m != j);
	}

	//Sort the eigen values (and vectors) in increase order
	for (int i = 0; i < n; i++) {
	int k = i;
	double p = realEigenvalues[i];
	for (int j = i + 1; j < n; j++) {
	if (realEigenvalues[j] > p) {
	k = j;
	p = realEigenvalues[j];
	}
	}
	if (k != i) {
	realEigenvalues[k] = realEigenvalues[i];
	realEigenvalues[i] = p;
	for (int j = 0; j < n; j++) {
	p = z[j][i];
	z[j][i] = z[j][k];
	z[j][k] = p;
	}
	}
	}

	// Determine the largest eigen value in absolute term.
	maxAbsoluteValue = 0;
	for (int i = 0; i < n; i++) {
	if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
	maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
	}
	}
	// Make null any eigen value too small to be significant
	if (maxAbsoluteValue!=0.0) {
	for (int i=0; i < n; i++) {
	if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
	realEigenvalues[i] = 0;
	}
	}
	}
	eigenvectors = new ArrayRealVector[n];
	final double[] tmp = new double[n];
	for (int i = 0; i < n; i++) {
	for (int j = 0; j < n; j++) {
	tmp[j] = z[j][i];
	}
	eigenvectors[i] = new ArrayRealVector(tmp);
	}
	}
	}