7
\$\begingroup\$

Following the 18.06 Linear algebra course, I was curious to reinvent the matrix class and basic functionality, like \$PA=LU \$ decomposition, Gauss elimination, finding inverse matrix etc.

Any comments about the code design, linear algebra and performance are welcome.

public class Matrix : ICloneable, IComparable<Matrix>
{
    private readonly double[,] _data;
    public int N => _data.GetUpperBound(0) + 1;
    public int M => _data.GetUpperBound(1) + 1;

    public Matrix(int n, bool diagonal = false)
    {
        _data = new double[n,n];
        if (!diagonal) return;
        for (int i = 0; i < n; i++)
        {
            _data[i, i] = 1.0;
        }
    }

    public Matrix(int n, int m, Random r = null)
    {
        _data = new double[n, m];
        if (r == null) return;
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < m; j++)
            {
                _data[i, j] = r.NextDouble();
            }
        }
    }

    public Matrix(double[,] data)
    {
        _data = data;
    }

    public ref double this[int row, int column] => ref _data[row, column];

    public static Matrix operator *(Matrix a, Matrix b)
    {
        if (a.M != b.N)
        {
            return null;
        }
        Matrix c = new Matrix(a.N, b.M);
        for (int i = 0; i < c.N; i++)
        {
            for (int j = 0; j < c.M; j++)
            {
                double s = 0.0;
                for (int m = 0; m < a.M; m++)
                {
                    s += a[i, m] * b[m, j];
                }
                c[i, j] = s;
            }
        }
        return c;
    }

    public static Matrix operator +(Matrix a, Matrix b)
    {
        if (a.M != b.M || a.N != b.N)
        {
            return null;
        }
        Matrix c = new Matrix(a.N, b.M);
        for (int i = 0; i < c.N; i++)
        {
            for (int j = 0; j < c.M; j++)
            {
                c[i, j] = a[i, j] + b[i, j];
            }
        }
        return c;
    }

    public void Transpose()
    {
        for (int i = 0; i < N; i++)
        {
            for (int j = i + 1; j < M; j++)
            {
                double tmp = _data[i, j];
                _data[i, j] = _data[j, i];
                _data[j, i] = tmp;
            }
        }
    }

    public void Transpose(out Matrix m)
    {
        m = new Matrix(N, M);
        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < M; j++)
            {
                m[i, j] = _data[j, i];
            }
        }
    }

    /// <summary>
    /// PA = LU factorization
    /// </summary>
    /// <param name="l">low-triangular matrix</param>
    /// <param name="p">permutation matrix</param>
    /// <param name="u">upper-triangular matrix</param>
    /// <returns></returns>
    public int PALU_factorization(out Matrix l, out Matrix p, out Matrix u)
    {
        return L_GaussEliminationForward(this, out l, out p, out u);
    }

    /// <summary>
    /// EPA = U factorization
    /// </summary>
    /// <param name="e">elimination matrix</param>
    /// <param name="p">permutation matrix</param>
    /// <param name="u">upper-triangular matrix</param>
    /// <returns></returns>
    public int EPAU_factorization(out Matrix e, out Matrix p, out Matrix u)
    {
        return E_GaussEliminationForward(this, out e, out p, out u);
    }

    /// <summary>
    /// Find out reverse matrix using Gauss-Jordan elimination
    /// </summary>
    /// <param name="reverseMatrix">reverse matrix to self</param>
    /// <returns>0 if success</returns>
    public int Reverse(out Matrix reverseMatrix)
    {
        if (N != M)
        {
            reverseMatrix = null;
            return -1;
        }
        int stdout = E_GaussEliminationForward(this, out var e, out var p, out var u);
        if (stdout != 0)
        {
            reverseMatrix = null;
            return stdout;
        }
        GaussEliminationBackward(u, e * p, out reverseMatrix);
        return 0;
    }

    /// <summary>
    /// Solve set of linear equations Ax=b
    /// </summary>
    /// <param name="b">right side matrix</param>
    /// <param name="x">matrix of variables</param>
    /// <returns>0 if success</returns>
    public int GaussElimination(Matrix b, out Matrix x)
    {
        if (N != b.N)
        {
            x = null;
            return -1;
        }
        int stdout = E_GaussEliminationForward(this, out var e, out var p, out var u);
        if (stdout != 0)
        {
            x = null;
            return stdout;
        }
        GaussEliminationBackward(u, e*p*b, out x);
        return 0;
    }

    /// <summary>
    /// Forward Gaussian Elimination to find E and P matrices from EPA = U equation
    /// </summary>
    /// <param name="a">coefficient matrix</param>
    /// <param name="e">eliminitaion matrix</param>
    /// <param name="p">permutation matrix</param>
    /// <param name="u">upper-triangular matrix</param>
    /// <returns>0 if success</returns>
    private static int E_GaussEliminationForward(Matrix a, out Matrix e, out Matrix p, out Matrix u)
    {
        e = new Matrix(a.N, true);
        p = new Matrix(a.N, true);
        u = (Matrix)a.Clone();
        for (int i = 0; i < a.N; i++)
        {
            if (Math.Abs(u[i, i]) < Double.Epsilon)
            {
                int iReverse = i;
                for (int j = i + 1; j < a.N; j++)
                {
                    if (Math.Abs(u[j, i]) > Double.Epsilon)
                    {
                        iReverse = j;
                        break;
                    }
                }

                if (iReverse == i)
                {
                    return -1;
                }
                e.ExchangeRows(iReverse, i, i);
                p.ExchangeRows(iReverse, i);
                u.ExchangeRows(iReverse, i);
            }
            Matrix eTmp = new Matrix(a.N, true);
            for (int j = i+1; j < a.N; j++)
            {
                double coeff = u[j, i] / u[i, i];
                eTmp[j, i] = -coeff;
                for (int k = i; k < a.M; k++)
                {
                    u[j, k] -= u[i, k] * coeff;
                }
            }
            e = eTmp*e;
        }
        return 0;
    }

    /// <summary>
    /// Forward Gaussian Elimination to find L and P matrices from PA = LU equation
    /// </summary>
    /// <param name="a">coefficient matrix (n*m)</param>
    /// <param name="l">low-triangular matrix (n*n)</param>
    /// <param name="p">permutation matrix (n*n)</param>
    /// <param name="u">upper-triangular matrix (n*m)</param>
    /// <returns>0 if success</returns>
    private static int L_GaussEliminationForward(Matrix a, out Matrix l, out Matrix p, out Matrix u)
    {
        l = new Matrix(a.N, true);
        p = new Matrix(a.N, true);
        u = (Matrix)a.Clone();
        for (int i = 0; i < a.N; i++)
        {
            if (Math.Abs(u[i, i]) < Double.Epsilon)
            {
                int iReverse = i;
                for (int j = i + 1; j < a.N; j++)
                {
                    if (Math.Abs(u[j, i]) > Double.Epsilon)
                    {
                        iReverse = j;
                        break;
                    }
                }

                if (iReverse == i)
                {
                    return -1;
                }
                l.ExchangeRows(iReverse, i, i);
                p.ExchangeRows(iReverse, i);
                u.ExchangeRows(iReverse, i);
            }
            for (int j = i + 1; j < a.N; j++)
            {
                double coeff = u[j, i] / u[i, i];
                l[j, i] = coeff;
                for (int k = i; k < a.M; k++)
                {
                    u[j, k] -= u[i, k] * coeff;
                }
            }
        }
        return 0;
    }

    /// <summary>
    /// Transforming augmented matrix [U|B] => [I|mB]
    /// </summary>
    /// <param name="u">upper-triangular matrix</param>
    /// <param name="c">right-side matrix</param>
    /// <param name="x">desired matrix</param>
    private void GaussEliminationBackward(Matrix u, Matrix c, out Matrix x)
    {
        x = (Matrix)c.Clone();
        for (int i = x.N - 1; i >= 0; i--)
        {
            for (int j = 0; j < x.M; j++)
            {
                x[i, j] /= u[i, i];
            }
            for (int j = 0; j < i; j++)
            {
                double coeff = u[j, i];
                for (int k = 0; k < x.M; k++)
                {
                    x[j, k] -= coeff * x[i, k];
                }
            }
        }
    }

    /// <summary>
    /// Exchange i, j rows of this.matrix until j-th column
    /// </summary>
    /// <param name="i">first row</param>
    /// <param name="j">second row</param>
    /// <param name="until">last column</param>
    public void ExchangeRows(int i, int j, int until = Int32.MaxValue)
    {
        until = Math.Min(M, until);
        for (int k = 0; k < until; k++)
        {
            double tmp = _data[i, k];
            _data[i, k] = _data[j, k];
            _data[j, k] = tmp;
        }
    }

    public override string ToString()
    {
        StringBuilder sb = new StringBuilder();
        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < M; j++)
            {
                sb.Append($"{_data[i, j]:0.000}\t");
            }
            sb.Append("\n");
        }
        return sb.ToString();
    }

    public int CompareTo(Matrix other)
    {
        if (N != other.N || M != other.M)
        {
            return -1;
        }

        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < M; j++)
            {
                if (Math.Abs(_data[i, j] - other[i, j]) > 0.0000000001)
                {
                    return -1;
                }
            }
        }

        return 0;
    }

    public object Clone()
    {
        return new Matrix((double[,])_data.Clone());
    }
}

Test

class Program
{
    static void Main(string[] args)
    {
        Matrix a = new Matrix(new double[,] {{1, 1, 1}, {2, 2, 5}, {4, 6, 8}});
        Console.WriteLine(a);
        a.Reverse(out var aRev);
        Console.WriteLine((a * aRev).CompareTo(aRev * a) == 0);
        Console.WriteLine((a * aRev).CompareTo(new Matrix(a.N, true)) == 0);
        a.PALU_factorization(out var l0, out var p0, out var u0);
        a.EPAU_factorization(out var e1, out var p1, out var u1);
        Console.WriteLine(p0.CompareTo(p1) == 0);
        Console.WriteLine(u0.CompareTo(u1) == 0);
        l0.Reverse(out var lRev);
        Console.WriteLine(lRev.CompareTo(e1) == 0);
        a.Transpose(out var aTra);
        aTra.Transpose();
        Console.WriteLine(a.CompareTo(aTra) == 0);
        Matrix b = new Matrix(new double[,] {{3}, {9}, {18}});
        a.GaussElimination(b, out var x);
        Console.WriteLine((a * x).CompareTo(b) == 0);
    }
}

Output

1.000   1.000   1.000
2.000   2.000   5.000
4.000   6.000   8.000

True
True
True
True
True
True
True
\$\endgroup\$
3
\$\begingroup\$

In the book “Introduction to Linear Algebra”, that comes along with this course, the author highlights several times the advantages of finding L (low-triangle matrix) instead of E (elimination matrix). It becomes pretty clear by inspecting the code.

private static int E_GaussEliminationForward(Matrix a, out Matrix e, out Matrix p, out Matrix u)
{
    // all necessary initialization
    e = new Matrix(a.N, true);
    for (int i = 0; i < a.N; i++)
    {
        // if pivot is 0, go and find an exchange
        Matrix eTmp = new Matrix(a.N, true);
        // elimination step
        e = eTmp*e; // TOO EXPENSIVE
    }
}

But I did not expect that the difference was going to be so huge. Just as an example, I implemented the Inverse function through finding matrix E and matrix L, and compare performance by BenchmarkDotNet with square random matrices (N = {5, 10, 20, 30, 40, 50}).

Barplot of means

public int ReverseE(out Matrix reverseMatrix)
{
    if (N != M)
    {
        reverseMatrix = null;
        return -1;
    }
    int stdout = E_GaussEliminationForward(this, out var e, out var p, out var u); // EPA = U
    if (stdout != 0)
    {
        reverseMatrix = null;
        return stdout;
    }
    UpperGaussEliminationBackward(u, e * p, out reverseMatrix); // Ux = EPI
    return 0;
}

public int ReverseL(out Matrix reverseMatrix)
{
    if (N != M)
    {
        reverseMatrix = null;
        return -1;
    }
    int stdout = L_GaussEliminationForward(this, out var l, out var p, out var u); // PA = LU
    if (stdout != 0)
    {
        reverseMatrix = null;
        return stdout;
    }
    // No should solve Ux = L^{-1}PI
    LowerGaussEliminationBackward(l, p, out var c); // Lc = PI
    UpperGaussEliminationBackward(u, c, out reverseMatrix); // Ux = c
    return 0;
}

Yes, we are doing additional LowerGaussEliminationBackward, but it does not influence a lot on performance. Even EPA=U decomposition is better to be done through L matrix.

This is the current state of the implementation.

\$\endgroup\$
2
\$\begingroup\$

API

I don't like that + and * operators return null on mismatched parameters: that's exactly the kind of thing that will obscure bugs and causes crashes 'down the line'. I'd have these throwing an exception explaining why the method failed (e.g. "Could not multiply matrices because first.M and second.N did not match"). This would crash violently as soon as the problem is detected, meaning you are close to the problem, and would tell you what the problem is. The same goes for Reverse.

I also don't understand why you have methods returning 0 on success, with -1 being a seemingly arbitrary alternative. If 'success' and 'failure' are the only conditions, then a boolean would do (which makes testing cleaner), or else use a dedicate return type or enum. I'm also not sure that stdout is a good name for the success flags. if (stdout != 0) is pretty cryptic.

It's great that you have inline documentation on some of your methods, but they should really document the 'return null on failure' conditions, and even N and M (and especially the random constructor) could do with documentation. I'd definitely want documentation on the Matrix(double[,] data) constructor, as it is presently unclear whether data is copied or not.

CompareTo

I don't like that CompareTo doesn't define an ordering: you are performing equality here, so implementing IEquatable<Matrix> would probably make more sense.

Clone

I'd suggest providing a Clone method that returns a copy of the matrix, and provide an explicit implementation for ICloneable.Clone. Copying a matrix requires effort or a cast at the moment, neither of which make programmers happy.

At the very least, it would remove code like u = (Matrix)a.Clone(); from your own methods, which just feels slightly unclean.

Transpose

These methods currently assume a square matrix, and will both crash on a non-square matrix. I'm not sure what the inplace version should do, but the other should presumably produce an M*N matrix.

Misc

  • In operator *, I wouldn't use m to count up to M: that's just asking for someone to make a mistake, and makes the code harder to read quickly.

  • Personally I'd write the 'random' constructor as 2 constructors, and I'd replace the 'diagonal' fielded constructor with a static Identity(int) method (or at least rename the diagonal to identity, which is more explicit): it's not clear what p = new Matrix(a.N, true) does at the moment.

  • Is there a particular reason you use GetUpperBound(i) + 1 rather than GetLength(i) to determine N and M?

  • I'm not qualified to comment on the quality of the numerical methods, but I'm a bit suspicious that you compare < Double.Epsilon and > Double.Epsilon, but never allow equality (of what is a perfectly valid value).

  • It seems odd to provide a public Reverse(out Matrix) method but no in-place equivalent given you provide both for Transpose. I also wonder why where is no operator - to go with operator +.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.