Click here to Skip to main content
11,710,943 members (81,812 online)
Click here to Skip to main content

Matrix And Polynominal Algebra

, 8 Apr 2014 CPOL 8.7K 372 43
Rate this:
Please Sign up or sign in to vote.
The article shows implementation of fundamental algorithms in matrix and polynominal algebra

Introduction

The fundamental algorithms used in matrix and polynominal algebra are to be presented in this article. All algorithms' implementations in polynominal are shown both in the fields of complex numbers. Algorithms' implementation in matrix (and vector) algebra are shown in the field of real numbers.

Background

The intermediate algebra definitions and operations knowledge will be required to follow implemenatation.

These are:

  • complex numbers operations;
  • basic opertaions at matrices (addition, subtraction, multiplication, scaling, transpose, trace, determinant);
  • matrix characteristic polynominal calculation using the Faddeev-Leverrier algorithm;
  • inverse of matrix calculation using the Faddeev-Leverrier algorithm;
  • polynominal derivative calculation algorithm;
  • the Horner's method of polynominal by binominal division and polynominal evaluation;
  • the Laguerre’s method for calculating polynominal roots;

The complex number structure and operations are provided with System.Numerics namespace.

Using the code

Convert Double to Complex method

The very simple methods that converts Double values to Complex structure.

public Complex Double2Complex(Double Value)
{
    Complex ValueComplex = new Complex(Value, 0);

    return ValueComplex;
}

 
public Complex[] Double2Complex(Double[] Values)
{
    Complex[] ValuesComplex = new Complex[Values.GetLength(0)];

    for (int m = 0; m < Values.GetLength(0); m++)
    {
        ValuesComplex[m] = new Complex(Values[m], 0);
    }

    return ValuesComplex;
} 
 

public Complex[,] Double2Complex(Double[,] Values)
{
    Complex[,] ValuesComplex = new Complex[Values.GetLength(0), Values.GetLength(1)];

    for (int m = 0; m < Values.GetLength(0); m++)
    {
        for (int n = 0; n < Values.GetLength(1); n++)
        {
            ValuesComplex[m, n] = new Complex(Values[m, n], 0);
        }
    }

    return ValuesComplex;
} 

Matrix operations

Identity matrix

The method returns square identity matrix of given dimension.

public Double[,] IdentityMatrix(int Dimension)
{
    if (Dimension > 0)
    {
        Double[,] I = new Double[Dimension, Dimension];
 
        for (int m = 0; m < Dimension; m++)
        {
            for (int n = 0; n < Dimension; n++)
            {
                if (m == n)
                {
                    I[m, n] = 1;
                }
                else
                {
                    I[m, n] = 0;
                }
            }
        }
 
        return I;
    }
    else
    {
        throw new Exception("Output identity matrix dimension must be positive.");
    }
}  

Matrix addition

The method returns a result of metrices addition. Both input matrices must have the same dimensions.

public Double[,] MatrixAddition(Double[,] MatrixA, Double[,] MatrixB)
{
    int M_A = MatrixA.GetLength(0);
    int N_A = MatrixA.GetLength(1);
    int M_B = MatrixB.GetLength(0);
    int N_B = MatrixB.GetLength(1);
 
    if (M_A == M_B && N_A == N_B)
    {
        Double[,] Sum = new Double[M_A, N_A];
 
        for (int m = 0; m < M_A; m++)
        {
            for (int n = 0; n < N_A; n++)
            {
                Sum[m, n] = MatrixA[m, n] + MatrixB[m, n];
            }
        }
 
        return Sum;
    }
    else
    {
        throw new Exception("Matrices dimensions must be equal.");
    }
}   

Matrix subtraction

The method returns a result of matrices subtraction. Both input matrices must have the same dimensions.

public Double[,] MatrixSubtraction(Double[,] MatrixA, Double[,] MatrixB)
{
    int M_A = MatrixA.GetLength(0);
    int N_A = MatrixA.GetLength(1);
    int M_B = MatrixB.GetLength(0);
    int N_B = MatrixB.GetLength(1);
 
    if (M_A == M_B && N_A == N_B)
    {
        Double[,] Sub = new Double[M_A, N_A];
 
        for (int m = 0; m < M_A; m++)
        {
            for (int n = 0; n < N_A; n++)
            {
                Sub[m, n] = MatrixA[m, n] - MatrixB[m, n];
            }
        }
 
        return Sub;
    }
    else
    {
        throw new Exception("Matrices dimensions must be equal.");
    }
}  

Matrix multiplication

The method returns a result of metrices multiplication. The number of columns in first multiplication component must be equal to the number of rows in second one.

public Double[,] MatrixMultiplication(Double[,] MatrixA, Double[,] MatrixB)
{
    int M_A = MatrixA.GetLength(0);
    int N_A = MatrixA.GetLength(1);
    int M_B = MatrixB.GetLength(0);
    int N_B = MatrixB.GetLength(1);
 
    if (N_A != M_B)
    {
        throw new Exception("Number of columns in A matrix must be equal to number of rows in B matrix.");
    }
    else
    {
        Double[,] _MatrixAB = new Double[M_A, N_B];
 
        for (int m_a = 0; m_a < M_A; m_a++)
        {
            for (int n_b = 0; n_b < N_B; n_b++)
            {
                _MatrixAB[m_a, n_b] = 0;
 
                for (int n_a = 0; n_a < N_A; n_a++)
                {
                    _MatrixAB[m_a, n_b] = _MatrixAB[m_a, n_b] + MatrixA[m_a, n_a] * MatrixB[n_a, n_b];
                }
            }
        }
 
        return _MatrixAB;
    }
}  

Matrix scaling

The method returns a result of matrix muliplication by constans coefficient.

public Double[,] MatrixScaling(Double Coefficient, Double[,] MatrixA)
{
    int M_A = MatrixA.GetLength(0);
    int N_A = MatrixA.GetLength(1);
 
    for (int m = 0; m < M_A; m++)
    {
        for (int n = 0; n < N_A; n++)
        {
            MatrixA[m, n] = Coefficient * MatrixA[m, n];
        }
    }
 
    return MatrixA;            
}  

Transpose of the matrix

The method returns a result of matrix transpose operation.

public Double[,] MatrixTranspose(Double[,] Matrix) 
{
    int M = Matrix.GetLength(0);
    int N = Matrix.GetLength(1);
 
    Double[,] _MatrixT = new Double[N, M];
 
    for (int m = 0; m < M; m++)
    {
        for (int n = 0; n < N; n++)
        {
            _MatrixT[n, m] = Matrix[m, n];
        }
    }
 
    return _MatrixT;
} 

Matrix trace

The method returns a matrix trace value - the sum of all diagonal elements. To calculate trace, input matrix must be square.

public Double MatrixTrace(Double[,] Matrix)
{
    int M_A = Matrix.GetLength(0);
    int N_A = Matrix.GetLength(1);
 
    if (M_A ==  N_A)
    {
        Double Trace = 0;
 
        for (int m = 0; m < M_A; m++)
        {
            Trace += Matrix[m, m];
        }
 
        return Trace;
    }
    else
    {
        throw new Exception("Matrix must be square.");
    }
}  

Matrix determinant

A recursive method using Laplace expansion and algebraic complement calculation.

public Double MatrixDeterminant(Double[,] Matrix)
{
    int M = Matrix.GetLength(0);
    int N = Matrix.GetLength(1);
 
    if (M == N)
    {
        Double Det = 0;
        
        if (M == 1)
        {
            Det = Matrix[0, 0];
        }
        else if (M == 2)
        {
            Det = Matrix[0, 0] * Matrix[1, 1] - Matrix[0, 1] * Matrix[1, 0];
        }
        else
        {
            Double[,] AlgebraicComplement = new Double[M - 1, M - 1];
 
            for (int m = 0; m < M; m++)
            {
                int a = 0;
 
                for (int k = 1; k < M; k++)
                {
                    int b = 0;
 
                    for (int l = 0; l < M; l++)
                    {
                        if (l != m)
                        {
                            AlgebraicComplement[a, b] = Matrix[k, l];
                            b++;
                        }
                    }
                    a++;
                }
                Det += Math.Pow(-1, m) * Matrix[0, m] * MatrixDeterminant(AlgebraicComplement);
            }
 
        }
 
        return Det;
    }
    else
    {
        throw new Exception("Matrix must be square.");
    }
}  

Matrix characteristic polynominal

The method calculates matrix characteristic polynominal coefficients using Faddeev - Leverrier algorithm. the method returns vector with elements corresponding to next polynominal coefficients starting with free term.

For polynominal with given description:

the vector of coefficients will be:

The same polynominal's coefficients structure will be used in all presented methods.

public Double[] CharacteristicPolynominalCoefficients(Double[,] Matrix)
{
    int M = Matrix.GetLength(0);
    int N = Matrix.GetLength(1);

    if (M == N)
    {
        Double[] Coeffs = new Double[M + 1];
        Double[] CoeffsSorted = new Double[M + 1];
        Double[,] B;
        Double LastCoeff;

        Coeffs[0] = 1;
        
        B = Matrix;
        LastCoeff = MatrixTrace(B);
        Coeffs[1] = -LastCoeff;

        for(int m = 2; m < M + 1; m++)
        {
            B = MatrixMultiplication(Matrix, MatrixSubtraction(B, MatrixScaling(LastCoeff, IdentityMatrix(B.GetLength(0)))));
            LastCoeff = MatrixTrace(B) / m;
            Coeffs[m] = -LastCoeff;
        }

        for (int m = 0; m < M + 1; m++)
        {
            CoeffsSorted[m] = Coeffs[M - m];
        }

        return CoeffsSorted;
    }
    else
    {
        throw new Exception("Input data matrix must be square.");
    }
}  

The inverse of matrix

The method calculates the inverse of matrix using Faddeev-Leverrier algorithm.

public Double[,] MatrixInverse(Double[,] Matrix)
{
    int M = Matrix.GetLength(0);
    int N = Matrix.GetLength(1);

    if (M == N)
    {
        if (MatrixDeterminant(Matrix) != 0)
        {
            Double[,] LastB;
            Double[,] NextB;
            Double[,] Inv = new Double[M, N];
            Double LastCoeff;
            Double NextCoeff;

            LastB = Matrix;
            LastCoeff = MatrixTrace(LastB);

            for (int m = 2; m < M; m++)
            {
                LastB = MatrixMultiplication(Matrix, MatrixSubtraction(LastB, MatrixScaling(LastCoeff, IdentityMatrix(LastB.GetLength(0)))));
                LastCoeff = MatrixTrace(LastB) / m;
            }

            NextB = MatrixMultiplication(Matrix, MatrixSubtraction(LastB, MatrixScaling(LastCoeff, IdentityMatrix(LastB.GetLength(0)))));
            NextCoeff = MatrixTrace(NextB) / M;

            Inv = MatrixScaling(1 / NextCoeff, MatrixSubtraction(LastB, MatrixScaling(LastCoeff, IdentityMatrix(LastB.GetLength(0)))));

            return Inv;
        }
        else
        {
            throw new Exception("The matrix is not invertible over commutative ring.");
        }
    }
    else
    {
        throw new Exception("Input data matrix must be square.");
    }
} 

Polynominal operations

Polynominal evaluation

The method calculates polynomianl evaluation for given value.

public Complex PolynominalEvaluation(Complex[] Coefficients, Complex Value)
{
    int M = Coefficients.GetLength(0);

    if (M > 1)
    {
        Complex[] CoeffsSorted = new Complex[M];
        Complex Last;

        for (int m = 0; m < M; m++)
        {
            CoeffsSorted[m] = Coefficients[M - m - 1];
        }

        Last = CoeffsSorted[0];

        for (int m = 1; m < M; m++)
        {
            Last = CoeffsSorted[m] + Last * Value;
        }

        return Last;
    }
    else
    {
        return Coefficients[0];
    }
} 

Polynominal addition

The method calculates the sum of two polynominals (the sum of coefficients).

public Complex[] PolynominalAddition(Complex[] Poly1, Complex[] Poly2)
{
    int M = Poly1.GetLength(0);
    int N = Poly2.GetLength(0);
    int K = Math.Max(M, N);
    Complex[] Poly1Ext = new Complex[K];
    Complex[] Poly2Ext = new Complex[K];
    Complex[] Add = new Complex[K];
    
    for (int m = 0; m < M; m++)
    {
        Poly1Ext[m] = Poly1[m];
    }

    for (int n = 0; n < N; n++)
    {
        Poly2Ext[n] = Poly2[n];
    }

    for (int k = 0; k < K; k++)
    {
        Add[k] = Poly1Ext[k] + Poly2Ext[k];
    }

    return Add;
} 

Polynominal scaling

The method calculates the polynominal multiplication by value.

public Complex[] PolynominalScaling(Complex[] Polynominal, Complex Value)
{
    int M = Polynominal.GetLength(0);

    for (int m = 0; m < M; m++)
    {
        Polynominal[m] = Value * Polynominal[m];
    }

    return Polynominal;
}  

Polynominal subtraction

Polynominal subtraction can be calculated using polynominal addition and polynominal scaling methods.

Complex[] Result = PolynominalAddition(Polynominal1, PolynominalScaling(Polynominal2, -1));  

Polynominal by binominal division

The method calculates a coefficients of new polynominal which is a result of polynominal and binominal division. The Horner's method is used for calculation. This method is not returning the remainder.

public Complex[] PolynominalByBinominalDivision(Complex[] Coefficients, Complex BinominalRoot)
{
    int M = Coefficients.GetLength(0);

    if (M > 1)
    {
        int N = M - 1;
        Complex[] Quotient = new Complex[N];
        Complex[] QuotientSorted = new Complex[N];
        Complex[] CoeffsSorted = new Complex[M];
        Complex Last;

        for (int m = 0; m < M; m++)
        {
            CoeffsSorted[m] = Coefficients[M - m - 1];
        }

        Last = CoeffsSorted[0];
        Quotient[0] = Last;

        for (int m = 1; m < N; m++)
        {
            Last = CoeffsSorted[m] + BinominalRoot * Last;
            Quotient[m] = Last;
        }

        Complex Remainder = CoeffsSorted[M - 1] + BinominalRoot * Last;

        for (int n = 0; n < N; n++)
        {
            QuotientSorted[n] = Quotient[N - n - 1];
        }

        return QuotientSorted;
    }
    else
    {
        throw new Exception("Given coefficients number is not corresponding to correct polynominal structure.");
    }
} 

Polynominal derivative

The method calculates a derivative of given polynominal.

public Complex[] PolynominalDerivative(Complex[] Coefficients)
{
    int Degree = Coefficients.GetLength(0) - 1;
    int NumCoefs = Coefficients.GetLength(0);
    Complex[] Derivative = new Complex[NumCoefs - 1];

    if (Degree > 0)
    {
        for (int i = 1; i < NumCoefs; i++)
        {
            Derivative[i - 1] = i * Coefficients[i];
        }
    }
    else
    {
        return new Complex[1] { 0 };
    }

    return Derivative;
} 

Polynominal roots

The method calculates polynominal roots with given accuracy and seed root (if unknown any real number can be given) using Laguerre’s algorithm. The method calculation assume that every polynominal root is a complex numer.

public Complex[] PolynominalRoots(Complex[] Coefficients, Complex Seed, Double Accuracy)
{
    int N = Coefficients.GetLength(0) - 1;
    Complex[] Roots = new Complex[N];
    Complex[] OrginalCoeffs = Coefficients;
    int degree = N;

    for (int n = 0; n < N; n++)
    {
        while (true)
        {
            Complex _tmp0 = PolynominalEvaluation(PolynominalDerivative(Coefficients), Seed);
            Complex _tmp1 = (degree - 1) * Complex.Pow(PolynominalEvaluation(PolynominalDerivative(Coefficients), Seed), 2);
            Complex _tmp2 = degree * PolynominalEvaluation(Coefficients, Seed) * PolynominalEvaluation(PolynominalDerivative(PolynominalDerivative(Coefficients)), Seed);
            Complex _tmp3 = Complex.Pow((degree - 1) * (_tmp1 - _tmp2), 0.5);
            Complex _tmp4 = MaximizeResultAbsolute(_tmp0, _tmp3);

            Seed = Seed - (degree * PolynominalEvaluation(Coefficients, Seed)) / _tmp4;
            Complex result = PolynominalEvaluation(Coefficients, Seed);

            if (Complex.Abs(result) < Math.Abs(Accuracy))
            {
                Roots[n] = Seed;
                Coefficients = PolynominalByBinominalDivision(Coefficients, Seed);
                degree--;
                Random _Random = new Random();
                Seed = -10 + (_Random.NextDouble() * 15);
                break;
            }
        }
    }

    return Roots;
} 

Maximize absolute value of operation result

The method calculates the result of operation (addition or subtraction) which gives greater absolute value. Used in Laguerre’s method.

public Complex MaximizeResultAbsolute(Complex Value1, Complex Value2)
{
    if (Complex.Abs(Value1 + Value2) > Complex.Abs(Value1 - Value2))
    {
        return Value1 + Value2;
    }
    else
    {
        return Value1 - Value2;
    }
} 

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Jakub Szymanowski
Systems Engineer Enspirion LLC (ENERGA Group)
Poland Poland
After graduating with a degree of Master of Science, Engineer in Automation and Robotics (at the Gdansk University of Technology) began his career as a Service Engineer of Ships' Systems. He spent four years in this profession and had a lot of places all over the world visited.

In 2012 he started to work for Polish energy company ENERGA (ENERGA Agregator LLC, ENERGA Innowacje LLC and Enspirion LLC) as Systems Engineer and Solutions Architect. Thereafter he participated in few large-scale research projects in the field of Information Technology and Systems Integration for Energy & Power Management and Demand Response applications as a Solution Architect, Business Architect or Project Manager.

His private interests include Signal and Image Processing, Adaptive Algorithms and Theory of Probability.

You may also be interested in...

Comments and Discussions

 
SuggestionMatrix And Polynominal Algebra Pin
Alexander Sharykin11-May-14 21:27
memberAlexander Sharykin11-May-14 21:27 
GeneralRe: Matrix And Polynominal Algebra Pin
Jakub Szymanowski11-May-14 21:47
professionalJakub Szymanowski11-May-14 21:47 
QuestionUnit Tests? Pin
Marc Clifton21-Apr-14 8:43
protectorMarc Clifton21-Apr-14 8:43 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Terms of Use | Mobile
Web03 | 2.8.150819.1 | Last Updated 8 Apr 2014
Article Copyright 2014 by Jakub Szymanowski
Everything else Copyright © CodeProject, 1999-2015
Layout: fixed | fluid