PowerIteration.cs

using System;
using System.Linq;
using Utilities.Extensions;

namespace Algorithms.LinearAlgebra.Eigenvalue
{
    /// <summary>
    ///     Power iteration method - eigenvalue numeric algorithm, based on recurrent relation:
    ///     Li+1 = (A * Li) / || A * Li ||, where Li - eigenvector approximation.
    /// </summary>
    public static class PowerIteration
    {
        /// <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>
        public static (double eigenvalue, double[] eigenvector) Dominant(
            double[,] source,
            double[] startVector,
            double error = 0.00001)
        {
            if (source.GetLength(0) != source.GetLength(1))
            {
                throw new ArgumentException("The source matrix is not square-shaped.");
            }

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

            double eigenNorm;
            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);

            var eigenvalue = 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>
        public static (double eigenvalue, double[] eigenvector) Dominant(double[,] source, double error = 0.00001) =>
            Dominant(source, new Random().NextVector(source.GetLength(1)), error);
    }
}