Improve Matrices Operations Speed (Part 2)

This new post will describe, how to invert a matrix and which are the best methods. As mentioned in “How to Solve” paragraph, the conventional Algorithm for the inversion between matrices is the LU decomposition. This method is quite expensive in terms of speed but I know there is a parallel version too. I have tried to study it, but I am sorry, I didn’t accept the challenge to rewrite LU decomposition in parallel version, it needed time that I didn’t have and I am only a developer and not mathematician.

Fortunately, when I studied Strassen for multiplication there was on the book an interesting conclusion: With this “fast” matrix multiplication, Strassen also obtained a surprising result for matrix inversion.

This very intriguing phrase bring me to the following formula:

R1 = Inverse(a11)
R2 = a21 x R1
R3 = R1 x a12
R4 = a21 x R3
R5 = R4 - a22
R6 = Inverse(R5)
c12 = R3 x R6
c21 = R6 x R2
R7 = R3 x c21
c11 = R1 - R7
c22 = -R6

So what I learned was that it is possible to obtain the matrix (multiple of 2^n) inversion from some “lighter” operations. I implemented this algorithm and compared it with the Gauss Jordan one.
Comparing those two algorithms, I obtained, using my notebook, these results:

Matrices 512x512 of Complex (double, double)
GJ; 24sec
Strassen: 19sec
20% time less

Matrices 1024x1024 of Complex (double, double)
GJ; 3min 02sec
Strassen: 2min 05sec
38% time less!

Used Invertion Methods

Gauss Jordan LU decomposition

public static Complex[,] GaussJordanInvertion(Complex[,] mat, BackgroundWorker bw = null)
{
    //Gestione notifica progressione            
    int iPerc = 0;
    int iMax = mat.GetLength(0);
    int iDone = 0;

    Complex[,] res;
    res = MatrixGeneric.Create<complex>(iMax);
    res = (Complex[,])mat.Clone();

    int[,] indxc;
    indxc = MatrixGeneric.Create<int>(iMax, true); 
    int[,] indxr;
    indxr = MatrixGeneric.Create<int>(iMax, true);
    int[,] ipiv;
    ipiv = MatrixGeneric.Create<int>(iMax, true);
    Complex big;
    Complex pivinv;
    Complex dum;
    int iRow = 1;
    int iCol = 1;

    for (int i = 1; i <= iMax; i++)
        ipiv[i, 1] = 0;
    for (int i = 1; i <= iMax; i++)
    {
        big = 0;
        for (int j = 1; j <= iMax; j++)
            if (ipiv[j, 1] != 1)
                for (int k = 1; k <= iMax; k++)
                {
                    if (ipiv[k, 1] == 0)
                    {
                        if (Complex.Abs(res[j, k]) >= Complex.Abs(big))
                        {
                            big = Complex.Abs(res[j, k]);
                            iRow = j;
                            iCol = k;
                        }
                    }
                    else if (ipiv[k, 1] > 1)
                    {
                        //Error "Singular Matrix -1"
                    }
                }
        ++(ipiv[iCol, 1]);
        if (iRow != iCol)
        {
            for (int l = 1; l <= iMax; l++)
                Swap(res[iRow, l], res[iCol, l]);
        }
        indxr[i, 1] = iRow;
        indxc[i, 1] = iCol;
        if (res[iCol, iCol] == 0)
        {
            //Error "Singular Matrix -2"
        }
        pivinv = 1 / res[iCol, iCol];
        res[iCol, iCol] = 1;
        for (int l = 1; l <= iMax; l++)
            res[iCol, l] = res[iCol, l] * pivinv;
        for (int ll = 1; ll <= iMax; ll++)
            if (ll != iCol)
            {
                dum = res[ll, iCol];
                res[ll, iCol] = 0;
                for (int l = 1; l <= iMax; l++)
                {
                    res[ll, l] = res[ll, l] - res[iCol, l] * dum;
                }
            }        
    }

    for (int l = iMax; l > 0; l--)
    {
        if (indxr[l, 1] != indxc[l, 1])
            for (int k = 1; k <= iMax; k++)
                Swap(res[k, indxr[l, 1]], res[k, indxc[l, 1]]);
    }

    return res;
}

Strassen

public static Complex[,] Inverse(Complex[,] a, int iThreshold = 128)
{
    int iOriginalSize = a.GetLength(0);
    int iFmtSize = GetSize(a.GetLength(0));
    Complex[,] A = a;

    int iSize = A.GetLength(0);
    int iHalfSize = iSize / 2;

    if (iSize <= iThreshold)    {
        return GaussJordanInvertion(a);
    }

    Complex[,] A11 = null;
    Complex[,] A12 = null;
    Complex[,] A21 = null;
    Complex[,] A22 = null;

    Complex[,] R1 = null;
    Complex[,] R2 = null;
    Complex[,] R3 = null;
    Complex[,] R4 = null;
    Complex[,] R5 = null;
    Complex[,] R6 = null;
    Complex[,] R7 = null;

    Complex[,] C11 = null;
    Complex[,] C12 = null;
    Complex[,] C21 = null;
    Complex[,] C22 = null;


    A11 = GetMatrixPortion(A, iHalfSize, 1, 1);
    R1 = Inverse(A11);

    A21 = GetMatrixPortion(A, iHalfSize, iHalfSize + 1, 1);
    R2 = Mul_Strassen(A21, R1, iThreshold);

    A12 = GetMatrixPortion(A, iHalfSize, 1, iHalfSize + 1);
    R3 = Mul_Strassen(R1, A12, iThreshold);

    R4 = Mul_Strassen(A21, R3, iThreshold);

    A22 = GetMatrixPortion(A, iHalfSize, iHalfSize + 1, iHalfSize + 1);
    R5 = Subtract(R4, A22);

    R6 = Inverse(R5);

    C12 = Mul_Strassen(R3, R6, iThreshold);
    C21 = Mul_Strassen(R6, R2, iThreshold);
    R7 = Mul_Strassen(R3, C21, iThreshold);

    C11 = Subtract(R1, R7);
    C22 = Multiply(R6, -1);

    Complex[,] c = MatrixGeneric.Create<complex>(iSize, iSize);
    for (int i = 1; i <= iHalfSize; i++)
    {
        for (int j = 1; j <= iHalfSize; j++)
        {
            c[i, j] = C11[i, j];
            c[i, j + iHalfSize] = C12[i, j];
            c[i + iHalfSize, j] = C21[i, j];
            c[i + iHalfSize, j + iHalfSize] = C22[i, j];
        }
    }

    A11 = null;
    A12 = null;
    A21 = null;
    A22 = null;

    R1 = null;
    R2 = null;
    R3 = null;
    R4 = null;
    R5 = null;
    R6 = null;
    R7 = null;

    C11 = null;
    C12 = null;
    C21 = null;
    C22 = null;

    GC.Collect();

    return c;
}

Remember that the input matrix shoud be a 2^n multiple, if not! complete it with 0 and 1 on the diagonal.