PowerIteration.cs

usingSystem;
usingSystem.Linq;
usingUtilities.Extensions;

namespaceAlgorithms.LinearAlgebra.Eigenvalue;

/// <summary>
/// Power iteration method - eigenvalue numeric algorithm, based on recurrent relation:
/// Li+1 = (A * Li) / || A * Li ||, where Li - eigenvector approximation.
/// </summary>
publicstaticclassPowerIteration
{
/// <summary>
/// Returns approximation of the dominant eigenvalue and eigenvector of <paramref name="source" /> matrix.
/// </summary>
/// <list type="bullet">
/// <item>
/// <description>The algorithm will not converge if the start vector is orthogonal to the eigenvector.</description>
/// </item>
/// <item>
/// <description>The <paramref name="source" /> matrix must be square-shaped.</description>
/// </item>
/// </list>
/// <param name="source">Source square-shaped matrix.</param>
/// <param name="startVector">Start vector.</param>
/// <param name="error">Accuracy of the result.</param>
/// <returns>Dominant eigenvalue and eigenvector pair.</returns>
/// <exception cref="ArgumentException">The <paramref name="source" /> matrix is not square-shaped.</exception>
/// <exception cref="ArgumentException">The length of the start vector doesn't equal the size of the source matrix.</exception>
publicstatic(doubleEigenvalue,double[]Eigenvector)Dominant(
double[,]source,
double[]startVector,
doubleerror=0.00001)
{
if(source.GetLength(0)!=source.GetLength(1))
{
thrownewArgumentException("The source matrix is not square-shaped.");
}

if(source.GetLength(0)!=startVector.Length)
{
thrownewArgumentException(
"The length of the start vector doesn't equal the size of the source matrix.");
}

doubleeigenNorm;
double[]previousEigenVector;
double[]currentEigenVector=startVector;

do
{
previousEigenVector=currentEigenVector;
currentEigenVector=source.Multiply(
previousEigenVector.ToColumnVector())
.ToRowVector();

eigenNorm=currentEigenVector.Magnitude();
currentEigenVector=currentEigenVector.Select(x =>x/eigenNorm).ToArray();
}
while(Math.Abs(currentEigenVector.Dot(previousEigenVector))<1.0-error);

vareigenvalue=source.Multiply(currentEigenVector.ToColumnVector()).ToRowVector().Magnitude();

return(eigenvalue,Eigenvector:currentEigenVector);
}

/// <summary>
/// Returns approximation of the dominant eigenvalue and eigenvector of <paramref name="source" /> matrix.
/// Random normalized vector is used as the start vector to decrease chance of orthogonality to the eigenvector.
/// </summary>
/// <list type="bullet">
/// <item>
/// <description>The algorithm will not converge if the start vector is orthogonal to the eigenvector.</description>
/// </item>
/// <item>
/// <description>The <paramref name="source" /> matrix should be square-shaped.</description>
/// </item>
/// </list>
/// <param name="source">Source square-shaped matrix.</param>
/// <param name="error">Accuracy of the result.</param>
/// <returns>Dominant eigenvalue and eigenvector pair.</returns>
/// <exception cref="ArgumentException">The <paramref name="source" /> matrix is not square-shaped.</exception>
/// <exception cref="ArgumentException">The length of the start vector doesn't equal the size of the source matrix.</exception>
publicstatic(doubleEigenvalue,double[]Eigenvector)Dominant(double[,]source,doubleerror=0.00001)=>
Dominant(source,newRandom().NextVector(source.GetLength(1)),error);
}