Matrix class in C#

Question

I've been learning C# during my free time in the past months; before that, I was mostly writing Java, so the transition hasn't been too hard, but I've never had my code reviewed or read by someone else.

• Is there anything that could be improved, or maybe simplified?
• Am I doing anything that would be considered bad practice or 'smelly' ? If so, what should I be doing instead?
• I very rarely write comments; should I use them more, and specifically here, are there any methods or sections that could use some comments?

The class DMatrix (D for double, as I intend to make one that uses floats instead, and another one for complex numbers), as you can guess, represents a matrix, or an m by n array of numbers. It started earlier today as rewrite of another implementation I had made in Java months ago; once I had finished translated that old Java code, I tried to include all the commonly used operations I could think of. I also tried to make the class easier and more intuitive to use, while keeping the code simple and -hopefully- readable, and trying to make use of all the C#-specific features I've learned about.

I've also included the Misc.Utils class, which mostly contains some static methods for generating random numbers in different ranges or types. I spend most of my coding time reinventing (or re-implementing) the wheel, or writing things from scratch, as I find it to be a good way to learn and understand most concepts, while also allowing me to implement and play with some interesting math or math-related topics.

DMatrix.cs

using System;
using Mathlib.Misc;

namespace Mathlib.Matrix {
    public class DMatrix {
        private double[,] data;

        public int NbRows {
            get;
            private set;
        }

        public int NbCols {
            get;
            private set;
        }

        public double this[int i, int j] {
            get {
                return data[i, j];
            }
            private set {
                data[i, j] = value;
            }
        }

        public static int Precision {
            get;
            set;
        }

        static DMatrix() {
            Precision = 2;
        }

        public DMatrix(int n) {
            if (n < 1) {
                throw new Exception($"Cannot create matrix," +
                                    $" N should be greater or equal to 1.\n\tN = {n}");
            }

            NbRows = NbCols = n;
            Fill(0);
        }

        public DMatrix(int m, int n) {
            if (m < 1) {
                throw new Exception($"Cannot create matrix," +
                                    $" M should be greater or equal to 1.\n\tM = {m}");
            }
            if (n < 1) {
                throw new Exception($"Cannot create matrix," +
                                    $" N should be greater or equal to 1.\n\tN = {n}");
            }

            NbRows = m;
            NbCols = n;

            data = new double[m, n];
            Fill(0);
        }

        public DMatrix(int m, int n, double x) {
            if (m < 1) {
                throw new Exception($"Cannot create matrix," +
                                    $" M should be greater or equal to 1.\n\tM = {m}");
            }
            if (n < 1) {
                throw new Exception($"Cannot create matrix," +
                                    $" N should be greater or equal to 1.\n\tN = {n}");
            }

            NbRows = m;
            NbCols = n;

            data = new double[m, n];
            Fill(x);
        }

        public DMatrix(double[,] data) {
            NbRows = data.GetLength(0);
            NbCols = data.GetLength(1);

            this.data = new double[NbRows, NbCols];

            for (int i = 0; i < NbRows; i++) {
                for (int j = 0; j < NbCols; j++) {
                    this.data[i, j] = data[i, j];
                }
            }
        }

        public DMatrix(DMatrix A) : this(A.data) { }

        public void Fill(double x) {
            for (int i = 0; i < NbRows; i++) {
                for (int j = 0; j < NbCols; j++) {
                    data[i, j] = x;
                }
            }
        }

        public void Rand() {
            for (int i = 0; i < NbRows; i++) {
                for (int j = 0; j < NbCols; j++) {
                    data[i, j] = Utils.Rand();
                }
            }
        }

        public void Rand(double max) {
            for (int i = 0; i < NbRows; i++) {
                for (int j = 0; j < NbCols; j++) {
                    data[i, j] = Utils.Rand(max);
                }
            }
        }

        public void Rand(double min, double max) {
            for (int i = 0; i < NbRows; i++) {
                for (int j = 0; j < NbCols; j++) {
                    data[i, j] = Utils.Rand(min, max);
                }
            }
        }

        public void SwapRows(int row1, int row2) {
            for (int i = 0; i < NbCols; i++) {
                double temp = data[row1, i];
                data[row1, i] = data[row2, i];
                data[row2, i] = temp;
            }
        }

        public static DMatrix operator -(DMatrix A) {
            DMatrix res = new DMatrix(A.NbRows, A.NbCols);

            for (int i = 0; i < A.NbRows; i++) {
                for (int j = 0; j < A.NbCols; j++) {
                    res.data[i, j] = -A.data[i, j];
                }
            }

            return res;
        }

        public static DMatrix operator ~(DMatrix A) {
            return A.Transpose();
        }

        public static DMatrix operator +(DMatrix A, DMatrix B) {
            if (A.NbRows != B.NbRows || A.NbCols != B.NbCols) {
                throw new Exception($"Cannot compute sum, matrix dimensions do not match:\n\t" +
                                    $"A : {A.NbRows}x{A.NbCols}\n\tB : {B.NbRows}x{B.NbCols}");
            }

            DMatrix res = new DMatrix(A.NbRows, A.NbCols);

            for (int i = 0; i < A.NbRows; i++) {
                for (int j = 0; j < A.NbCols; j++) {
                    res.data[i, j] = A.data[i, j] + B.data[i, j];
                }
            }

            return res;
        }

        public static DMatrix operator -(DMatrix A, DMatrix B) {
            return A + (-B);
        }

        public static DMatrix operator *(DMatrix A, DMatrix B) {
            if (A.NbCols != B.NbRows) {
                throw new Exception($"Cannot compute product, matrix dimensions do not match:\n\t" +
                                    $"A : {A.NbRows}x{A.NbCols}\n\tB : {B.NbRows}x{B.NbCols}");
            }

            DMatrix res = new DMatrix(A.NbRows, B.NbCols);

            for (int i = 0; i < res.NbRows; i++) {
                for (int j = 0; j < res.NbCols; j++) {
                    for (int k = 0; k < A.NbCols; k++) {
                        res.data[i, j] += A.data[i, k] * B.data[k, j];
                    }
                }
            }

            return res;
        }

        public DMatrix HadamardProduct(DMatrix A) {
            if (NbCols != A.NbCols || NbRows != A.NbRows) {
                throw new Exception($"Cannot compute Hadamard product," +
                                    $" dimensions do not match:\n\tA : {NbRows}x{NbCols}\n\t" +
                                    $"B : {A.NbRows}x{A.NbCols}");
            }

            DMatrix res = new DMatrix(NbRows, NbCols);

            for (int i = 0; i < res.NbRows; i++) {
                for (int j = 0; j < res.NbCols; j++) {
                    res.data[i, j] = data[i, j] * A.data[i, j];
                }
            }

            return res;
        }

        public static DMatrix operator +(DMatrix A, double x) {
            DMatrix res = new DMatrix(A.NbRows, A.NbCols);

            for (int i = 0; i < A.NbRows; i++) {
                for (int j = 0; j < A.NbCols; j++) {
                    res.data[i, j] = A.data[i, j] + x;
                }
            }

            return res;
        }

        public static DMatrix operator -(DMatrix A, double x) {
            return A + (-x);
        }

        public static DMatrix operator *(DMatrix A, double x) {
            DMatrix res = new DMatrix(A.NbRows, A.NbCols);

            for (int i = 0; i < A.NbRows; i++) {
                for (int j = 0; j < A.NbCols; j++) {
                    res.data[i, j] = A.data[i, j] * x;
                }
            }

            return res;
        }

        public double Trace() {
            if (NbCols != NbRows) {
                throw new Exception($"Cannot compute trace, matrix is not square:\n\t" +
                                    $"A : {NbRows}x{NbCols}");
            }

            double res = 0;

            for (int i = 0; i < NbRows; i++) {
                res += data[i, i];
            }

            return res;
        }

        public double Det_recursive() {
            if (NbCols != NbRows) {
                throw new Exception($"Cannot compute determinant, matrix isn't square:\n\t" +
                                    $"A : {NbRows}x{NbCols}");
            }

            double sum = 0;

            if (NbRows == 1) {
                return data[0, 0];
            }
            else if (NbRows == 2) {
                return data[0, 0] * data[1, 1] - data[0, 1] * data[1, 0];
            }

            for (int N = 0; N < NbRows; N++) {
                DMatrix minor = new DMatrix(NbRows - 1, NbRows - 1);

                for (int i = 0; i < NbRows - 1; i++) {
                    for (int j = 0; j < NbRows; j++) {
                        if (j < N) {
                            minor.data[i, j] = data[i + 1, j];
                        }
                        else if (j > N) {
                            minor.data[i, j - 1] = data[i + 1, j];
                        }
                    }
                }

                int sgn = (N % 2 == 0) ? 1 : -1;
                sum += sgn * data[0, N] * minor.Det_recursive();
            }

            return sum;
        }

        public double Norm_frobenius() {
            double res = 0;

            if (NbRows == 1 || NbCols == 1) {
                for (int i = 0; i < NbRows; i++) {
                    for (int j = 0; j < NbCols; j++) {
                        res += data[i, j] * data[i, j];
                    }
                }
            }
            else {
                res = (~this * this).Trace();
            }

            return Math.Sqrt(res);
        }

        public static DMatrix Id(int n) {
            DMatrix res = new DMatrix(n);

            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    if (i == j) {
                        res.data[i, j] = 1;
                    }
                    else {
                        res.data[i, j] = 0;
                    }
                }
            }

            return res;
        }

        public DMatrix Transpose() {
            DMatrix res = new DMatrix(NbRows, NbCols);

            for (int i = 0; i < NbRows; i++) {
                for (int j = 0; j < NbCols; j++) {
                    data[i, j] = data[j, i];
                }
            }

            return res;
        }

        public DMatrix Row(int row) {
            DMatrix A = new DMatrix(1, NbCols);

            for (int i = 0; i < NbCols; i++) {
                A.data[0, i] = data[row, i];
            }

            return A;
        }

        public DMatrix Col(int col) {
            DMatrix A = new DMatrix(NbRows, 1);

            for (int i = 0; i < NbRows; i++) {
                A.data[i, 0] = data[i, col];
            }

            return A;
        }

        public DMatrix Diag() {
            if (NbRows != NbCols) {
                throw new Exception($"Cannot compute diagonal, matrix isn't square:\n\t" +
                                    $"A : {NbRows}x{NbCols}");
            }

            DMatrix res = new DMatrix(NbRows, NbCols);

            for (int i = 0; i < NbRows; i++) {
                res.data[i, i] = data[i, i];
            }

            return res;
        }

        public DMatrix Map(OneVarFunction func) {
            DMatrix res = new DMatrix(NbRows, NbCols);

            for (int i = 0; i < res.NbRows; i++) {
                for (int j = 0; j < res.NbCols; j++) {
                    res.data[i, j] = func(data[i, j]);
                }
            }

            return res;
        }

        public double[,] ToArray() {
            double[,] arr = new double[NbRows, NbCols];

            for (int i = 0; i < NbRows; i++) {
                for (int j = 0; j < NbCols; j++) {
                    arr[i, j] = data[i, j];
                }
            }

            return arr;
        }

        public string ToCsv() {
            string str = "";

            for (int i = 0; i < NbRows; i++) {
                for (int j = 0; j < NbCols; j++) {
                    if (j != NbCols - 1) {
                        str += $"{data[i, j]},";
                    }
                    else {
                        str += $"{data[i, j]}\n";
                    }
                }
            }

            return $"{str}\n";
        }

        public override string ToString() {
            string str = "";

            for (int i = 0; i < NbRows; i++) {
                str += "[ ";
                for (int j = 0; j < NbCols; j++) {
                    string decimals = new string('0', Precision);
                    string nb = data[i, j].ToString($"+0.{decimals};-0.{decimals};0")
                        .Replace('+', ' ');

                    if (j != NbCols - 1) {
                        str += nb + "  ";
                    }
                    else {
                        str += nb + " ]\n";
                    }
                }
            }

            return str + '\n';
        }
    }
}

Utils.cs

using System;

namespace Mathlib.Misc {
    public static class Utils {
        private static int seed;
        private static Random rng;

        public static void Seed(int s) {
            seed = s;
            rng = new Random(s);
        }

        public static int Seed() {
            return seed;
        }

        static Utils() {
            Seed((int) DateTime.Now.Ticks);

            rng = new Random(seed);
        }

        public static double Rand() {
            return rng.NextDouble();
        }

        public static int Rand(int max) {
            return (int) (rng.NextDouble() * max);
        }

        public static float Rand(float max) {
            return (float) (rng.NextDouble() * max);
        }

        public static double Rand(double max) {
            return rng.NextDouble() * max;
        }

        public static int Rand(int min, int max) {
            return (int) (rng.NextDouble() * (max - min) + min);
        }

        public static float Rand(float min, float max) {
            return (float) (rng.NextDouble() * (max - min) + min);
        }

        public static double Rand(double min, double max) {
            return rng.NextDouble() * (max - min) + min;
        }

        public static float Lerp(float val, float x1, float x2, float y1, float y2) {
            if (val == x1) {
                return y1;
            }

            if (val == x2) {
                return y2;
            }

            return (y2 - y1) / (x2 - x1) * (val - x1) + y1;
        }

        public static double Lerp(double val, double x1, double x2, double y1, double y2) {
            if (val == x1) {
                return y1;
            }

            if (val == x2) {
                return y2;
            }

            return (y2 - y1) / (x2 - x1) * (val - x1) + y1;
        }
    }
}
```

Answer 1

I like that your braces and indentation is beautifully consistent. But, this is C#, and your Java is showing ;-) - Kudos for the mostly-consistent PascalCase type and member names, but it's the K&R same-line { opening brace that clashes with the typical Allman next-line { opening brace standard most people would expect of C# code:

public DMatrix(int n)
{
    if (n < 1)
    {
        throw new Exception($"Cannot create matrix," +
                            $" N should be greater or equal to 1.\n\tN = {n}");
    }

    NbRows = NbCols = n;
    Fill(0);
}

I would write auto-properties on a single line:

public int NbRows { get; private set; }

public int NbCols { get; private set; }

Or rather, I'd make these get-only:

public int Rows => _rows;

public int Columns => _columns;

..and have the corresponding private fields. Now maybe I'm old-school, but I like that underscore for private fields: it removes the need to qualify private fields with this when you'd like a local variable by that name. But, something feels wrong about needing to store this kind of metadata in private fields... let's keep reading..

public static int Precision {
    get;
    set;
}

Uh-oh. Assuming it's used elsewhere and not just in the ToString implementation (why then?), it's almost guaranteed you don't want a static modifier here: a static property belongs to the type, not the object/instance - if I have two DMatrix objects and I say matrix1.Precision = 4 and then two lines later I say matrix2.Precision = 1, it's perfectly reasonable to expect matrix1.Precision to not start claiming its value is now 1... right?

public DMatrix(DMatrix A) : this(A.data) { }

Hm, I was about to write something about constructor chaining.. there's not enough of that:

public DMatrix() : this(1) { }

public DMatrix(int size) : this(size, size)

public DMatrix(int width, int height) : this(width, height, default) { }

public DMatrix(int width, int height, int fillValue) : this(new double[width, height](), fillValue) { }

public DMatrix(double[,] matrix, int fillValue = default)
{
    /* one place */
}

I'd make the encapsulated data array readonly (that makes the array reference read-only, not the array elements; the compiler will prevent assigning to the array outside a constructor).

    private readonly double[,] _data;

And then, the Rows and Cols don't need a setter anymore:

public int Rows { get; }
public int Columns { get; }

They can only be assigned in a constructor.

Careful with operator overloads - it's not clear how the rather obscure unary one's complement operator turns into an implicit shorthand for applying the ~ operator to each element ...transposing:

public static DMatrix operator ~(DMatrix A) {
    return A.Transpose();
}

There's already a Transpose method for that.

I see you throw Exception - avoid that: you want to throw a meaningful exception type, such that you know what the problem is just with the type name. When throwing in a guard clause that's validating the given arguments, you want to throw an ArgumentException for example.

if (NbCols != NbRows) {
    throw new Exception($"Cannot compute trace, matrix is not square:\n\t" +
                        $"A : {NbRows}x{NbCols}");
}

This particular condition should be handled in the constructor IMO (where NbCols/Columns and NbRows/Rows are assigned), and then here a System.Diagnostics.Debug.Assert call should suffice to cover this rather unlikely scenario - it's the job of the constructor(s) to ensure the object is in a valid state.

This one however:

if (m < 1) {
    throw new Exception($"Cannot create matrix," +
                        $" M should be greater or equal to 1.\n\tM = {m}");
}

...worries me. If NbRows = NbCols = n, and then data = new double[m, n];, it means our parameterization is 1-based, but the matrix itself is 0-based... so we ask for 1x1 and actually get a 2x2 where [0,0] is shoved under the proverbial carpet!

A note here:

public double Det_recursive() {

The common guidelines would scream for PascalCase here, and I would suggest that the recursive nature of the implementation is ..an implementation detail, and so I would have gone with Determinant for a name.

I like how you antagonized + and - operators, and I know you really like the K&R braces, but this:

public static DMatrix operator -(DMatrix A, double x) {
    return A + (-x);
}

could be:

public static DMatrix operator -(DMatrix A, double x) => A + (-x);

Look ma, no hands braces!

---

There's definitely other things to say, notably about Utils, its name and static nature, and how this impacts coupling, and how this coupling impacts testability... but I'll let other reviewers cover other aspects of this code.

Answer 2

A bug: Transpose does data[i, j] = data[j, i];, but it returns a new matrix, it's not supposed to change the matrix in-place, and that doesn't actually work (going out of range for a non-square matrix and losing data for a square one).

Manual loops to copy between arrays: did you know you can use Array.Copy on multi-dimensional arrays? The indexes work as if the rows are laid end-to-end.

Det_recursive: cofactor expansion is a well known example of an algorithm with O(n!) time complexity. It's fine for a tiny matrix yes, but a factorial quickly gets out of hand. Naming it Determinant is more dangrous if you leave it like this, offering no warning that it's going to be done with the slow algorithm, but fixing the efficiency is much nicer than putting a warning in the name. It can be done in O(n³) time via LU decomposition, which is something that's also useful for various other purposes.

Unnecessary zero filling: the constructors and the Id function wastes a bunch of time and code writing zeroes into a new array full of zeroes. Id can be simplified by just writing the ones to the diagonal, and the constructor without fill parameter could simply not fill.

Answer 3

New line character differs depending on system. Can't rely on \n. Let me introduce you to your new friend, Environment.NewLine

Concatenating strings in a loop is very inefficient. It's recommended to use StringBuilder instead

Code reuse. There are a couple of spots where duplicated code can be extracted away.

• Iterating an array and setting values. This can be extracted into a method than accepts a function (T currentValue, int x, int y) -> T newValue and sets each cell to the value returned by the function. Please see SetValues and its usages in the code below.
• Creating a new matrix. This has been abstracted into several methods in my code (Clone(), CreateNewMatrixWithSameSize() and CreateNewMatrix()).

Generics. You mention that you want to later implement matrixes for different kind of numbers. A good idea would be to use a generic abstract class Matrix<T>. I've done some of the work in the code below, but let me summarize some of the changes I needed to make:

• Change all mentions of double to T
• Extract all code that operates with doubles or instantiates matrixes to it's own abstract methods.
• Create DoubleMatrix class inheriting from Matrix<T> and implement the abstract methods we created in the previous step.

• Move other type specific logic (Like Rand() and Id()) to DoubleMatrix class

  Matrix.cs

public abstract class Matrix<T>
{
    protected Matrix(int n) : this(n, n)
    {
    }

    protected Matrix(int m, int n) : this(m, n, default(T))
    {
    }

    protected Matrix(int m, int n, T x)
    {
        if (m < 1)
        {
            throw new Exception($"Cannot create matrix," +
                                $" M should be greater or equal to 1.\n\tM = {m}");
        }

        if (n < 1)
        {
            throw new Exception($"Cannot create matrix," +
                                $" N should be greater or equal to 1.\n\tN = {n}");
        }

        NbRows = m;
        NbCols = n;

        data = new T[m, n];
        Fill(x);
    }

    protected Matrix(Matrix<T> A) : this(A.data)
    {
    }

    protected Matrix(T[,] source)
    {
        NbRows = source.GetLength(0);
        NbCols = source.GetLength(1);

        this.data = new T[NbRows, NbCols];

        this.SetValues((i, j) => source[i, j]);
    }

    private readonly T[,] data;

    public int NbRows { get; }
    public int NbCols { get; }

    public T this[int i, int j]
    {
        get => data[i, j];
        private set => data[i, j] = value;
    }

    public static Matrix<T> operator -(Matrix<T> a) => a.Negative();
    public static Matrix<T> operator ~(Matrix<T> a) => a.Transpose();
    public static Matrix<T> operator +(Matrix<T> a, Matrix<T> B) => a.Add(B);
    public static Matrix<T> operator -(Matrix<T> a, Matrix<T> B) => a + (-B);
    public static Matrix<T> operator *(Matrix<T> a, Matrix<T> B) => a.Multiply(B);
    public static Matrix<T> operator +(Matrix<T> a, T x) => a.Add(x);
    public static Matrix<T> operator -(Matrix<T> a, T x) => a.Subtract(x);
    public static Matrix<T> operator *(Matrix<T> a, T x) => a.MultiplyBy(x);

    public void Fill(T x) => SetValues(x);

    public void SwapRows(int row1, int row2)
    {
        for (int i = 0; i < NbCols; i++)
        {
            T temp = this[row1, i];
            this[row1, i] = this[row2, i];
            this[row2, i] = temp;
        }
    }

    public Matrix<T> HadamardProduct(Matrix<T> A)
    {
        if (NbCols != A.NbCols || NbRows != A.NbRows)
        {
            throw new Exception($"Cannot compute Hadamard product," +
                                $" dimensions do not match:\n\tA : {NbRows}x{NbCols}\n\t" +
                                $"B : {A.NbRows}x{A.NbCols}");
        }

        return CreateSameSizeMatrix()
            .SetValues((i, j) => Multiply(this[i, j], A[i, j]));
    }

    public T Trace()
    {
        if (NbCols != NbRows)
        {
            throw new Exception($"Cannot compute trace, matrix is not square:\n\t" +
                                $"A : {NbRows}x{NbCols}");
        }

        T res = default(T);

        for (int i = 0; i < NbRows; i++)
        {
            res = Add(res, this[i, i]);
        }

        return res;
    }

    public T Determinant()
    {
        if (NbCols != NbRows)
        {
            throw new Exception($"Cannot compute determinant, matrix isn't square:\n\t" +
                                $"A : {NbRows}x{NbCols}");
        }

        T sum = default(T);

        if (NbRows == 1)
        {
            return this[0, 0];
        }
        else if (NbRows == 2)
        {
            T member1 = Multiply(this[0, 0], this[1, 1]);
            T member2 = Multiply(this[0, 1], this[1, 0]);

            return Subtract(member1, member2);
        }

        for (int N = 0; N < NbRows; N++)
        {
            Matrix<T> minor = CreateMatrix(NbRows - 1, NbRows - 1);

            for (int i = 0; i < NbRows - 1; i++)
            {
                for (int j = 0; j < NbRows; j++)
                {
                    if (j < N)
                    {
                        minor[i, j] = this[i + 1, j];
                    }
                    else if (j > N)
                    {
                        minor[i, j - 1] = this[i + 1, j];
                    }
                }
            }

            T quotient = (N % 2 == 0) ? this[0, N] : Negative(this[0, N]);
            T product = Multiply(minor.Determinant(), quotient);
            sum = Add(sum, product);
        }

        return sum;
    }

    public T Norm_frobenius()
    {
        T res = default(T);

        if (NbRows == 1 || NbCols == 1)
        {
            for (int i = 0; i < NbRows; i++)
            {
                for (int j = 0; j < NbCols; j++)
                {
                    T product = Multiply(this[i, j], this[i, j]);
                    res = Add(res, product);
                }
            }
        }
        else
        {
            res = (~this * this).Trace();
        }

        return Sqrt(res);
    }

    public Matrix<T> Transpose() => Clone().SetValues((_, i, j) => this[j, i]);

    public Matrix<T> Row(int row)
    {
        Matrix<T> A = CreateMatrix(1, NbCols);

        for (int i = 0; i < NbCols; i++)
        {
            A[0, i] = this[row, i];
        }

        return A;
    }

    public Matrix<T> Column(int col)
    {
        Matrix<T> A = CreateMatrix(NbRows, 1);

        for (int i = 0; i < NbRows; i++)
        {
            A[i, 0] = data[i, col];
        }

        return A;
    }

    public Matrix<T> Diagonal()
    {
        if (NbRows != NbCols)
        {
            throw new Exception($"Cannot compute diagonal, matrix isn't square:\n\t" +
                                $"A : {NbRows}x{NbCols}");
        }

        Matrix<T> res = CreateMatrix(NbRows, NbCols);

        for (int i = 0; i < NbRows; i++)
        {
            res[i, i] = this[i, i];
        }

        return res;
    }

    public Matrix<T> Map(Func<T, T> func) =>
        CreateMatrix(NbRows, NbCols)
            .SetValues(func);

    public T[,] ToArray()
    {
        T[,] arr = new T[NbRows, NbCols];

        for (int i = 0; i < NbRows; i++)
        {
            for (int j = 0; j < NbCols; j++)
            {
                arr[i, j] = this[i, j];
            }
        }

        return arr;
    }

    public string ToCsv()
    {
        var sb = new StringBuilder();

        for (int i = 0; i < NbRows; i++)
        {
            for (int j = 0; j < NbCols; j++)
            {
                if (j != NbCols - 1)
                {
                    sb.Append($"{this[i, j]},");
                }
                else
                {
                    sb.Append($"{this[i, j]}{Environment.NewLine}");
                }
            }
        }

        return sb.ToString();
    }

    public override string ToString()
    {
        var sb = new StringBuilder();

        for (int i = 0; i < NbRows; i++)
        {
            sb.Append("[ ");
            for (int j = 0; j < NbCols; j++)
            {

                string nb = AsString(this[i, j]);

                if (j != NbCols - 1)
                {
                    sb.Append($"{nb} ");
                }
                else
                {
                    sb.Append($"{nb} ]{Environment.NewLine}");
                }
            }
        }

        return sb.ToString();
    }

    protected abstract Matrix<T> Clone();
    protected abstract Matrix<T> CreateMatrix(int m, int n);
    private Matrix<T> CreateSameSizeMatrix() => CreateMatrix(NbRows, NbCols);
    protected abstract T Negative(T x);
    protected abstract T Add(T x, T y);
    protected abstract T Multiply(T x, T y);
    protected abstract T Sqrt(T x);
    protected virtual string AsString(T x) => x.ToString();
    private T Subtract(T x, T y) => Add(x, Negative(y));

    private Matrix<T> Multiply(Matrix<T> B)
    {
        if (NbCols != B.NbRows)
        {
            throw new Exception($"Cannot compute product, matrix dimensions do not match:\n\t" +
                                $"A : {NbRows}x{NbCols}\n\tB : {B.NbRows}x{B.NbCols}");
        }

        Matrix<T> res = CreateMatrix(NbRows, B.NbCols);

        res.SetValues((current, i, j) =>
        {
            var sum = default(T);
            for (int k = 0; k < NbCols; k++)
            {
                sum = Add(sum, Multiply(this[i, k], B[k, j]));
            }

            return sum;
        });

        return res;
    }

    private Matrix<T> Negative() => Clone().SetValues(Negative);

    private Matrix<T> Add(Matrix<T> B)
    {
        if (NbRows != B.NbRows || NbCols != B.NbCols)
        {
            throw new Exception($"Cannot compute sum, matrix dimensions do not match:\n\t" +
                                $"A : {NbRows}x{NbCols}\n\tB : {B.NbRows}x{B.NbCols}");
        }

        return Clone()
            .SetValues((current, i, j) => Add(current, B[i, j]));
    }

    private Matrix<T> Add(T x) =>
        CreateSameSizeMatrix()
            .SetValues((i, j) => Add(this[i, j], x));

    private Matrix<T> Subtract(T x) => Add(Negative(x));

    private Matrix<T> MultiplyBy(T x) =>
        CreateSameSizeMatrix()
            .SetValues((i, j) => Multiply(this[i, j], x));

    protected Matrix<T> SetValues(Func<T, int, int, T> newValueSelector)
    {
        for (int i = 0; i < NbRows; i++)
        {
            for (int j = 0; j < NbCols; j++)
            {
                this[i, j] = newValueSelector(this[i, j], i, j);
            }
        }

        return this;
    }

    protected Matrix<T> SetValues(T newValue) =>
        SetValues((_, __, ___) => newValue);

    protected Matrix<T> SetValues(Func<T, T> newValueSelector) =>
        SetValues((current, _, __) => newValueSelector(current));

    protected Matrix<T> SetValues(Func<int, int, T> newValueSelector) =>
        SetValues((_, i, j) => newValueSelector(i, j));
}

DoubleMatrix.cs

public class DoubleMatrix : Matrix<double>
{
    public DoubleMatrix(int n) : base(n) {}
    public DoubleMatrix(int m, int n) : base(m, n) {}
    public DoubleMatrix(int m, int n, double x) : base(m, n, x) {}
    public DoubleMatrix(double[,] source) : base(source) {}
    public DoubleMatrix(Matrix<double> A) : base(A) {}

    public static DoubleMatrix Identity(int n) => 
        (DoubleMatrix) new DoubleMatrix(n).SetValues((i, j) => i == j ? 1 : 0);

    protected override Matrix<double> Clone() => new DoubleMatrix(this);
    protected override Matrix<double> CreateMatrix(int m, int n) => new DoubleMatrix(m, n);
    protected override double Negative(double x) => -x;
    protected override double Add(double x, double y) => x + y;
    protected override double Multiply(double x, double y) => x * y;
    protected override double Sqrt(double x) => Math.Sqrt(x);
}

Answer 4

Next step: Make it Expression<>-oriented.

Making your own math tools is tons of fun and practical!

Since you're asking about potential improvements, my suggestion is to move toward Expression<>-oriented coding next.

First, you define:

public abstract partial class Expression<T>
{
    public T Evaluate()
    {
        return this.Internal_Evaluate();
    }
    protected abstract T Internal_Evaluate();
}

You can ignore the .Evaluate()/.Internal_Evaluate() distinction for now, though I'd suggest that you include it as it may make your life easier later.

Anyway, then you can define stuff like constants

public partial class ConstantExpression<T>
    :        Expression<T>
{
    protected T ConstantValue { get; private set; }

    //  Empty constructor that no external code should ever use:
    protected ConstantExpression() { }

    //  Primary factory-style constructor.
    //  If you add overloads, try to make them call this one,
    //  such that this is the only method that ever includes
    //  "new ConstantExpression<>()" anywhere in your code.
    public static ConstantExpression<T> New(
                T constantValue
        )
    {
        var toReturn = new ConstantExpression<T>();

        toReturn.ConstantValue = constantValue;

        System.Threading.Thread.MemmoryBarrier();  // Just always include this until you have a reason not to.

        return toReturn;
    }

    protected override T Internal_Evaluate()
    {
        return this.ConstantValue;
    }
}

and addition

public partial class AdditionExpression
    :        Expression<double>
{
    protected Expression<double> Argument0Expression { get; private set; }
    protected Expression<double> Argument1Expression { get; private set; }

    //  Empty constructor that no external code should ever use:
    protected AdditionExpression() { }

    //  Primary factory-style constructor.
    //  If you add overloads, try to make them call this one,
    //  such that this is the only method that ever includes
    //  "new AdditionExpression()" anywhere in your code.
    public static AdditionExpression New(
                Expression<double> argument0Expression
            ,   Expression<double> argument1Expression
        )
    {
        if (argument0Expression == null || argument1Expression == null)
        {
            throw new Exception();  // replace with your preferred debug-tracing style
        }

        var toReturn = new AdditionExpression();

        toReturn.Argument0Expression = argument0Expression;
        toReturn.Argument1Expression = argument1Expression;

        System.Threading.Thread.MemmoryBarrier();  // Just always include this until you have a reason not to.

        return toReturn;
    }

    protected override double Internal_Evaluate()
    {
        var argument0 = this.Argument0Expression.Evaluate();
        var argument1 = this.Argument1Expression.Evaluate();

        return argument0 + argument1;
    }
}

with usability helpers like

partial class Expression<T>
{
    public static implicit operator Expression<T>(
                T constantValue
        )
    {
        return ConstantExpression<T>.New(constantValue);
    }

    public static Expression<double> operator +(
                Expression<double> addend0Expression
            ,   Expression<double> addend1Expression
        )
    {
        return AdditionExpression.New(
                   addend0Expression
               ,   addend1Expression
            );
    }
}

Then now that you've got the basic outline for Expression<>'s, you can rewrite your matrix code:

public partial class MatrixExpression
    :        Expression<double[,]>
{
    protected Expression<double>[,] MatrixElementsExpressions { get; private set; }

    protected MatrixExpression() { }
    public static MatrixExpression New(
                Expression<double[,]> matrixElementsExpressions
        )
    {
        var toReturn = new MatrixExpression();

        toReturn.MatrixElementsExpressions = matrixElementsExpressions;

        System.Threading.Thread.MemoryBarrier();

        return toReturn;
    }

    protected override double[,] Internal_Evaluate()
    {
        var matrixElementsExpressions = this.MatrixElementsExpressions;

        var length_0 = matrixElementsExpressions.GetLength(0);
        var length_1 = matrixElementsExpressions.GetLength(1);

        var toReturn = new double[length_0, length_1];

        for (long i_0 = 0; i_0 < length_0; ++i_0)
        {
            for (long i_1 = 0; i_1 < length_1; ++i_1)
            {
                toReturn[i_0, i_1] = matrixElementsExpressions[i_0, i_1].Evaluate();
            }
        }

        return toReturn;
    }
}

Then, it might be tempting to add, say, a .Transpose() method to MatrixExpresion – but don't!

Instead:

public static Expression<double[,]> Transpose(
            this Expression<double[,]> matrixExpression
    )
{
    if (matrixExpression == null)
    {
        throw new Exception();    //  Replace with your preferred error-handling system.
    }

    var toReturn = TransposedMatrixExpression.New(
                matrixExpression
        );

    return toReturn;
}

public partial class TransposedMatrixExpression
    :        Expression<double[,]>
{
    protected Expression<double[,]> MatrixExpression { get; private set; }

    protected TransposedMatrixExpression() { }
    public static TransposedMatrixExpression New(
                Expression<double[,]> matrixExpression
        )
    {
        var toReturn = new TransposedMatrixExpression();

        toReturn.MatrixExpression = matrixExpression;

        System.Threading.Thread.MemoryBarrier();

        return toReturn;
    }

    protected override double[,] Internal_Evaluate()
    {
        var matrixExpression = this.MatrixExpression;

        var matrix = matrixExpression.Evaluate();

        var length_0 = matrix.GetLength(0);
        var length_1 = matrix.GetLength(1);

        var toReturn = new double[length_1, length_0];

        for (long i_0 = 0; i_0 < length_0; ++i_0)
        {
            for (long i_1 = 0; i_1 < length_1; ++i_1)
            {
                toReturn[i_1, i_0] = matrix[i_0, i_1];
            }
        }

        return toReturn;
    }
}

In general, keep Expression<>-definitions slim. New operations shouldn't be additional methods within other classes, but rather each get its own Expression<>, e.g. as we effectively added a .Transpose() method via the class TransposedMatrixExpression above.

Note: Classes and methods are the same thing.

This may be confusing, so I'm quote-boxing it out: you can ignore this point if it doesn't make sense.

C# methods and C# classes are logically equivalent, if we ignore some variation in presentation and implied implementation details. To better understand this, you might look into how anonymous C# methods get their own C# class in the runtime.

Once you understand their equivalence, it'll help frame why we define methods as classes, e.g. as with .Transpose() above.

---

Tips

1. Put operator definitions, e.g. +, into a partial class Expression<>{ } block, as we did for +(Expression<double>, Expression<double>) above.

2. My coding style may look verbose. I've left a lot of room for things that I suspect most people will want to add as they develop a project like this. I'd advise against trying to make it shorter for a long time, until you get several steps beyond this.

3. It may feel weird to have ConstantExpression<T>'s getting trivially .Evaluate()'d to the T .ConstantValue that they wrap. You may feel inclined to try to reduce overhead by, for example, defining a variant of AdditionExpression that works on a T and an Expression<T>, rather than two Expression<T>'s, to help reduce unnecessary method calls. If you feel strongly about trying this, then it can be a good learning experience – but, it's a mistake.

---

Next steps.

Obviously, there's a lot to play with here. That's a lot of fun!

Then you can also do stuff like:

1. Creating graphical interfaces for Expression<>'s.

   • I originally did this with WPF. I'd suggest coding it purely in C#; ignore the XML interface.

2. Add in symbolic logic.

3. Extend the logic beyond math into general programming structures.

4. Add in calculus, differential equations, etc..

5. Hoist the whole thing onto a custom evaluation engine.

   • At first, stick with just having an Expression<>-tree doing depth-first .Evaluate()-ing, as above.