# Matrix And Polynominal Algebra

0 1,027

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 eval‎uation;

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. Collapse | Copy Code

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. Collapse | Copy Code

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.”);

}

}

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

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. Collapse | Copy Code

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. Collapse | Copy Code

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. Collapse | Copy Code

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. Collapse | Copy Code

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. Collapse | Copy Code

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. Collapse | Copy Code

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. Collapse | Copy Code

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 = 1;

B = Matrix;

LastCoeff = MatrixTrace(B);

Coeffs = -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. Collapse | Copy Code

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 eval‎uation

The method calculates polynomianl eval‎uation for given value. Collapse | Copy Code

public Complex PolynominalEval‎uation(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;

for (int m = 1; m < M; m++)

{

Last = CoeffsSorted[m] + Last * Value;

}

return Last;

}

else

{

return Coefficients;

}

}

The method calculates the sum of two polynominals (the sum of coefficients). Collapse | Copy Code

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];

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++)

{

}

}

Polynominal scaling

The method calculates the polynominal multiplication by value. Collapse | Copy Code

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. Collapse | Copy Code

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. Collapse | Copy Code

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;

Quotient = 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. Collapse | Copy Code

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 { 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. Collapse | Copy Code

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 = PolynominalEval‎uation(PolynominalDerivative(Coefficients), Seed);

Complex _tmp1 = (degree – 1) * Complex.Pow(PolynominalEval‎uation(PolynominalDerivative(Coefficients), Seed), 2);

Complex _tmp2 = degree * PolynominalEval‎uation(Coefficients, Seed) * PolynominalEval‎uation(PolynominalDerivative(PolynominalDerivative(Coefficients)), Seed);

Complex _tmp3 = Complex.Pow((degree – 1) * (_tmp1 – _tmp2), 0.5);

Complex _tmp4 = MaximizeResultAbsolute(_tmp0, _tmp3);

Seed = Seed – (degree * PolynominalEval‎uation(Coefficients, Seed)) / _tmp4;

Complex result = PolynominalEval‎uation(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. Collapse | Copy Code

public Complex MaximizeResultAbsolute(Complex Value1, Complex Value2)

{

if (Complex.Abs(Value1 + Value2) > Complex.Abs(Value1 – Value2))

{

return Value1 + Value2;

}

else

{

return Value1 – Value2;

}

}