A Second Wind for Square Matrix Algebra Thanks to New Algorithms

Introduction: Rethinking Square Matrix Functions

Matrix analysis serves as a powerful tool for generalizing, unifying, systematizing, and classifying a wide variety of mathematical objects. It finds applications in algorithm analysis, digital signal processing, economics, biology, and many other fields. At the same time, it remains a core analytical framework for theoretical research.

Among the fundamental functions associated with matrices are the determinant and the permanent. However, these classical functions come with major challenges. To this day, no general recurrence relations (formulas that express terms of a sequence in terms of their predecessors and indices) have been established for either the determinant or the permanent. Moreover, no polynomial-time algorithm is known for computing the permanent of a square matrix, and a general representation of these functions in terms of partition polynomials remains an open problem.

Can we construct a function on square matrices whose structure allows it to be computed using linear recurrence relations, and which also satisfies the two fundamental properties that define the determinant?

Recent developments suggest that this is indeed possible. A new class of matrix functions has been introduced that not only meets these criteria but also opens up new directions for both theoretical and applied research. Preliminary results show strong connections with known mathematical sequences, support efficient algorithmic computation, and have potential applications in fields such as symbolic computation, combinatorics, systems modeling, and cryptography.

These functions may offer a new foundation for algorithmic mathematics — one that leverages the power of modern computing to handle tasks previously considered too complex. For the IT sector in particular, this could mean faster algorithms for matrix evaluation, enhanced tools for data modeling, and novel approaches to problems involving large-scale discrete structures.

Reinventing Matrix Algebra: A Fresh Take on an Old Structure

Matrix functions like the determinant and the permanent are deeply rooted in the structure of permutations. These classical tools evaluate square matrices by summing over permutations of their entries — but what if we could go further?

A new class of matrix functions, introduced in 2025, builds upon this foundation by using the concept of transversals — sets of matrix entries that contain exactly one element from each row and column. This approach allows for a richer structure and deeper combinatorial interpretation.

To define this new function — the par-function — several conditions were imposed:

each term of the function’s multilinear polynomial is the product of entries from a transversal of the matrix;
every matrix entry influences the value of the function;
half of the terms in the polynomial appear with a positive sign, while the other half appear with a negative sign.

Since both the determinant and the permanent of an 𝑛×𝑛 square matrix are defined based on the set of all permutations of the first 𝑛 positive integers, it is natural to expect that the new multilinear matrix function will exhibit certain connections with the determinant and the permanent.

Let 𝐴 be a square matrix of order 𝑛:

$A_n = \left( \begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{array} \right)$

We define the par-function of the matrix 𝐴 as the multilinear polynomial:

$par(A_n) = \sum_{r=1}^{n} \sum_{(i_1,i_2,\ldots,i_n)_r \in CS_n} (-1)^{n-r} a_{1i_1} a_{2i_2} \cdots a_{ni_n}$

where 𝐶𝑆𝑛 denotes a certain subset of the symmetric group 𝑆𝑛, the set of all permutations of the first 𝑛 positive integers.

There exists a bijection (one-to-one correspondence) between the set of all permutations of the first 𝑛 positive integers and the set of ordered partitions of the number 𝑛.

To illustrate, we present the entries of a matrix of order 4 that will be used to construct the terms of the par-function according to the ordered partitions of the number 4:

$Custom Web App Development$

Analogously to the definition of the permanent, we introduce the par+ function, in which all terms are taken with a positive sign:

$\text{par}^{+}(A_n) = \sum_{r=1}^{n} \sum_{(i_1, i_2, \ldots, i_n)_r \in CS_n} a_{1i_1} a_{2i_2} \cdots a_{ni_n}$

The following recurrence relations will then hold:

$\text{par}(A_n) = \sum_{i=1}^{n} \prod_{j=0}^{i-1} (-1)^{i-1} a_{n-j,\;n-i+j+1} \, \text{par}(A_{n-i})$

$\text{par}^{+}(A_n) = \sum_{i=1}^{n} \prod_{j=0}^{i-1} a_{n-j,\;n-i+j+1} \; \text{par}^{+}(A_{n-i})$

The functions described above were first introduced in 2025 by the author of this article, a programmer-mathematician, in the scientific paper “par-Functions of Square Matrices.”

From Theory to Code: Computer Algorithms for Par and Par+

To investigate the properties of the new functions of square matrices, it will be most advantageous to develop computer algorithms. We will provide some examples in C#, as it is particularly well-suited for parallel computation, which is crucial for the analysis.

1. Algorithm for computing the par and par+ functions:

public static class PFuncAlg
{ 
             /// <summary>
             /// Calculates par and par+ functions for a square integer matrix.
             /// </summary>
    public static (int pFunc, int pPlusFunc) CalculatePFunctionsForMatrix(int[,] matrix)
    {
        var size = (byte)matrix.GetLength(0);
        var partitions = CalculatePartitions(size);
        return CalculatePFunctions(matrix, partitions, size);
    }

             /// <summary>
             /// Generates all ordered partitions of a given number.
             /// </summary>
    public static List<byte[]> CalculatePartitions(byte size)
    {
        var result = new List<byte[]>();
        var currentPartition = new List<byte>();

        void GeneratePartitions(byte remaining)
        {
            currentPartition.Add(0);
            var lastIndex = currentPartition.Count - 1;

            for (byte i = 1; i <= remaining; i++)
            {
                currentPartition[lastIndex] = i;

                if (remaining == i)
                {
                    result.Add([.. currentPartition]);
                }
                else
                {
                    GeneratePartitions((byte)(remaining - i));
                }
            }
            currentPartition.RemoveAt(lastIndex);
        }
        GeneratePartitions(size);
        return result;
    }
    
             /// <summary>
             /// Evaluates par and par+ using the given partitions.
             /// </summary>
    public static (int pFunc, int pPlusFunc) CalculatePFunctions(int[,] matrix, List<byte[]> partitions, byte size)
    {
        var pFunc = 0;
        var pPlusFunc = 0;

        foreach (var partition in partitions)
        {
            var part = GetPartitionValue(matrix, partition);
            if (part != 0)
            {
                // Alternating sign based on parity of (size - partition length)
                if ((size - partition.Length) % 2 == 0)
                {
                    pFunc += part;
                }
                else
                {
                    pFunc -= part;
                }
                pPlusFunc += part;
            }
        }
        return (pFunc, pPlusFunc);
    }
    
             /// <summary>
             /// Computes product of matrix entries defined by a partition.
             /// </summary>
    private static int GetPartitionValue(int[,] matrix, byte[] partition)
    {
        var res = 1;
        var x = 0;
        var y = -1;

        foreach (var p in partition)
        {
            var oldY = y;
            y += p;
            for (var i = y; i > oldY; i--)
            {
                var e = matrix[x, i];
                if (e == 0)
                {
                    return 0;
                }
                res *= e;
                x++;
            }
        }
        return res;
    }
}

2. Algorithm for computing “weight” and min/max of par and par+ for binary matrices:

public static class PFuncAnalyzerAlg
{
    // Parallel options for max degree of parallelism based on available processors
    private static readonly ParallelOptions ParallelOptionsMax = new() { MaxDegreeOfParallelism = Environment.ProcessorCount };

    // Lock object to ensure thread-safety during critical sections
    private static readonly object Lock = new();
    
             /// <summary>
             /// Calculates the gravity matrices (par-function and par+ function) for a given matrix size.
             /// </summary>
    public static (int[,], int[,]) GetGravityMatrix(byte size)
    {
        var partitions = PFuncAlg.CalculatePartitions(size);
        var primaryCombinations = GetAllPossibleCombinations(size, 1);

        var pRes = GetMatrix([], size);
        var pPlusRes = GetMatrix([], size);

        Parallel.ForEach(primaryCombinations, ParallelOptionsMax, combination =>
        {
            var matrix = GetMatrix(combination, size);
            var (pFunc, pPlusFunc) = PFuncAlg.CalculatePFunctions(matrix, partitions, size);
            pRes[combination[0].x, combination[0].y] = pFunc;
            pPlusRes[combination[0].x, combination[0].y] = pPlusFunc;
        });
        return (pRes, pPlusRes);
    }
    
             /// <summary>
             /// Calculates the minimum and maximum values for par-function and par+ function based on a specified number of zero elements in the matrix.
             /// </summary>
    public static (int, int, int, int) GetMinMax(byte size, byte amountOfZeroElements)
    {
        var partitions = PFuncAlg.CalculatePartitions(size);
        var primaryCombinations = GetAllPossibleCombinations(size, amountOfZeroElements);

        var pMin = int.MaxValue;
        var pMax = int.MinValue;
        var pPlusMin = int.MaxValue;
        var pPlusMax = int.MinValue;

        Parallel.ForEach(primaryCombinations, ParallelOptionsMax, pCombination =>
        {
            var matrix = GetMatrix(pCombination, size);
            var (pFunc, pPlusFunc) = PFuncAlg.CalculatePFunctions(matrix, partitions, size);

            Monitor.Enter(Lock);
            if (pFunc > pMax) pMax = pFunc;
            if (pFunc < pMin) pMin = pFunc;
            if (pPlusFunc > pPlusMax) pPlusMax = pPlusFunc;
            if (pPlusFunc < pPlusMin) pPlusMin = pPlusFunc;
            Monitor.Exit(Lock);
        });
        return (pMin, pMax, pPlusMin, pPlusMax);
    }
    
             /// <summary>
             /// Creates a matrix with elements set to zero at given positions.
             /// </summary>
    private static int[,] GetMatrix((byte x, byte y)[] elements, byte size)
    {
        var res = new int[size, size];
        for (byte i = 0; i < size; i++)
        {
            for (byte j = 0; j < size; j++)
            {
                res[i, j] = 1;
            }
        }
        foreach (var (x, y) in elements)
        {
            res[x, y] = 0;
        }
        return res;
    }
    
             /// <summary>
             /// Generates all possible combinations of positions in the matrix where a specified number of elements will be set to zero.
             /// </summary>
    private static IEnumerable<(byte x, byte y)[]> GetAllPossibleCombinations(byte size, byte elementsAmount)
    {
        if (elementsAmount == 0)
            return [];

        var totalPositions = size * size;
        var allPositions = new (byte x, byte y)[totalPositions];

        var index = 0;
        for (byte x = 0; x < size; x++)
        {
            for (byte y = 0; y < size; y++)
            {
                allPositions[index++] = (x, y);
            }
        }
        var maxIndexBase = totalPositions - elementsAmount;
        var combination = new (byte x, byte y)[elementsAmount];

        IEnumerable<(byte x, byte y)[]> GenerateCombinations(int startIndex, int depth)
        {
            var nextDepth = depth + 1;
            var maxIndex = maxIndexBase + nextDepth;
            if (maxIndex == totalPositions)
            {
                for (var i = startIndex; i < maxIndex; i++)
                {
                    combination[depth] = allPositions[i];
                    var res = new (byte x, byte y)[elementsAmount];
                    Array.Copy(combination, res, elementsAmount);
                    yield return res;
                }
            }
            else
            {
                for (var i = startIndex; i < maxIndex; i++)
                {
                    combination[depth] = allPositions[i];
                    foreach (var sub in GenerateCombinations(i + 1, nextDepth))
                    {
                        yield return sub;
                    }
                }
            }
        }
        return GenerateCombinations(0, 0);
    }
}

Expanding Horizons: Future Research and Cross-Disciplinary Applications

The introduction of the par-function has opened the door to classifying a wide array of recurrence relations that appear across mathematical disciplines — a feat not previously almost achievable. Using computer-aided analysis, the par-function has been shown to connect with numerous known integer sequences and has even led to the discovery of several new ones.

Moreover, the par-function has practical applications in solving systems of ordinary differential equations, which model phenomena in physics, economics, biology, and other domains. By leveraging the structure of fundamental solution systems, it offers new ways to analyze and solve these equations.

Ongoing research aims to further explore these directions. The author is currently preparing scientific papers that will detail these and other emerging applications of the par and par+ functions. Future work will also explore their relevance to other areas of science, where matrix-based structures and transformations play a key role.

Let’s Summarize Everything Together

Mathematics offers a universal language for describing the processes that govern our universe. Among its most powerful tools are matrices, differential equations, and recurrence relations. This article has explored how one can, through the dual lenses of programming and mathematics, define and analyze entirely new matrix functions.

These functions not only inherit key properties of well-known constructs like the determinant and the permanent but also include characteristics with promising applications. Such functions have been systematically constructed and studied using algorithmic methods. Leveraging computational power, this research enables the exploration of mathematical ideas at a scale and speed that would have been unattainable using traditional approaches—potentially accelerating progress in several fields of modern science.

Innovation Through Science at Diatom Enterprises

At Diatom Enterprises, we are proud to support and encourage this kind of scientific exploration by our developers. Innovation is not just a byproduct of our engineering process — it is a core part of our culture. Whether through pioneering algorithmic research or creating high-performance business solutions, we believe that pushing the boundaries of knowledge leads to exceptional products. Our commitment extends beyond writing clean, efficient code: we strive to cultivate a space where deep theoretical thinking and practical development come together, advancing both technology and science.

Table of contents

math tutoring programming