/*
* LargestEigenvalueFinder.java
*
* Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard
*
* This file is part of BEAST.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* BEAST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*/
package dr.math.matrixAlgebra;
import dr.math.iterations.IterativeProcess;
/**
* Object used to find the largest eigen value and the corresponding
* eigen vector of a matrix by successive approximations.
*
* @author Didier H. Besset
*/
public class LargestEigenvalueFinder extends IterativeProcess
{
/**
* Eigenvalue
*/
private double eigenvalue;
/**
* Eigenvector
*/
private Vector eigenvector;
/**
* Eigenvector of transposed matrix
*/
private Vector transposedEigenvector;
/**
* Matrix.
*/
private Matrix matrix;
/**
* Constructor method.
* @param prec double
* @param a MatrixAlgebra.Matrix
*/
public LargestEigenvalueFinder ( double prec, Matrix a)
{
this(a);
this.setDesiredPrecision ( prec);
}
/**
* Constructor method.
* @param a MatrixAlgebra.Matrix
*/
public LargestEigenvalueFinder ( Matrix a)
{
matrix = a;
eigenvalue = Double.NaN;
}
/**
* Returns the eigen value found by the receiver.
* @return double
*/
public double eigenvalue ( )
{
return eigenvalue;
}
/**
* Returns the normalized eigen vector found by the receiver.
* @return MatrixAlgebra.Vector
*/
public Vector eigenvector ( )
{
return eigenvector.product( 1.0 / eigenvector.norm());
}
/**
* Iterate matrix product in eigenvalue information.
*/
public double evaluateIteration()
{
double oldEigenvalue = eigenvalue;
transposedEigenvector =
transposedEigenvector.secureProduct( matrix);
transposedEigenvector = transposedEigenvector.product( 1.0
/ transposedEigenvector.components[0]);
eigenvector = matrix.secureProduct( eigenvector);
eigenvalue = eigenvector.components[0];
eigenvector = eigenvector.product( 1.0 / eigenvalue);
return Double.isNaN( oldEigenvalue)
? 10 * getDesiredPrecision()
: Math.abs( eigenvalue - oldEigenvalue);
}
/**
* Set result to undefined.
*/
public void initializeIterations()
{
eigenvalue = Double.NaN;
int n = matrix.columns();
double [] eigenvectorComponents = new double[ n];
for ( int i = 0; i < n; i++) { eigenvectorComponents [i] = 1.0;}
eigenvector = new Vector( eigenvectorComponents);
n = matrix.rows();
eigenvectorComponents = new double[ n];
for ( int i = 0; i < n; i++) { eigenvectorComponents [i] = 1.0;}
transposedEigenvector = new Vector( eigenvectorComponents);
}
/**
* Returns a finder to find the next largest eigen value of the receiver's matrix.
* @return MatrixAlgebra.LargestEigenvalueFinder
*/
public LargestEigenvalueFinder nextLargestEigenvalueFinder ( )
{
double norm = 1.0 / eigenvector.secureProduct(
transposedEigenvector);
Vector v1 = eigenvector.product( norm);
return new LargestEigenvalueFinder( getDesiredPrecision(),
matrix.secureProduct(SymmetricMatrix.identityMatrix(
v1.dimension()).secureSubtract(v1.tensorProduct(
transposedEigenvector))));
}
/**
* Returns a string representation of the receiver.
* @return java.lang.String
*/
public String toString()
{
StringBuffer sb = new StringBuffer();
sb.append( eigenvalue);
sb.append(" (");
sb.append( eigenvector.toString());
sb.append(')');
return sb.toString();
}
}