Immutable Matrix

Question

I'm writing implementations of some numerical methods to solve linear equations systems, those implementations use the following Matrix class. I'm trying to get this class immutable, due to that the constructor and some methods create a copies of the backing two-dimesional array. The class have some basic functionality to operate NxM Matrix e.g. overloaded operators for sum, difference and product. The methods Equals(object) and HashCode() were stubbed by Resharper. I'd would like to read any suggestions and improvements.

 public class Matrix
{              
    public int M { get; private set; }
    public int N { get; private set; }        
    private readonly double[,] _values;               

    public Matrix(double[,] values)
    {            
        M = values.GetLength(0);
        N = values.GetLength(1);
        _values = (double[,])values.Clone();            
    }

    public double this[int row, int column] {
        get
        {
            ValidateRowIndex(row);
            ValidateColumnIndex(column);
            return _values[row, column];
        }
    }                 

    public IEnumerable<double> Column(int column)
    {
        ValidateColumnIndex(column);

        for (int i = 0; i <= _values.GetUpperBound(0); ++i)
        {
            yield return _values[i, column];
        }
    }

    public IEnumerable<double> Row(int row)
    {
        ValidateRowIndex(row);

        for (int j = 0; j <= _values.GetUpperBound(1); ++j)
        {
            yield return _values[row, j];
        }
    }

    /// <summary>
    /// Method that will add a row times a factor to other row
    /// in the system.
    /// </summary>
    /// <param name="row1">Row to be added</param>
    /// <param name="row2">Row that will be added to</param>
    /// <param name="factor">Times the row1 will be multiplied before be added</param>
    /// <returns>Result matrix</returns>
    public Matrix AddRow(int row1, int row2, double factor)
    {
        ValidateRowIndex(row1);
        ValidateRowIndex(row2);
        if (factor == 0)
        {
            throw new ArgumentException("factor cannot be zero");
        }

        var result = new double[M, N];
        var row = Row(row1).ToArray();

        for (int i = 0; i < M; i++)
        {
            for (int j = 0; j < N; j++)
            {
                result[i, j] = (i != row2) ?
                                  this[i, j] :
                                  this[i, j] + (row[j] * factor);
            }                
        }

        return result;
    }

    public static Matrix operator +(Matrix a, Matrix b)
    {
        ValidateSameDimensions(a, b);

        return Operation(a, b, (n1, n2) => n1 + n2);
    }

    public static Matrix operator -(Matrix a, Matrix b)
    {
        ValidateSameDimensions(a, b);

        return Operation(a, b, (n1, n2) => n1 - n2);
    }

    public static Matrix operator *(Matrix a, Matrix b)
    {
        if (a.M != b.N || a.N != b.M)
        {
            throw new ArgumentException("Matrix are not product compatible");
        }

        var result = new double[a.M, b.N];

        for (int i = 0; i < a.M; i++)
        {
            for (int j = 0; j < b.N; j++)
            {
                for (int k = 0; k < b.M; k++)
                {
                    result[i, j] = result[i, j] + (a[i, k] * b[k, j]);
                }
            }
        }

        return result;
    }

    public static Matrix operator *(double scalar, Matrix m)
    {
        var result = new double[m.M, m.N];
        for (int i = 0; i < m.M; i++)
        {
            for (int j = 0; j < m.N; j++)
            {
                result[i, j] = scalar * m[i, j];
            }
        }

        return result;
    }

    public static Matrix operator *(Matrix m, double scalar)
    {
        return scalar * m;
    }

    public static implicit operator Matrix(double[,] values)
    {
        return new Matrix(values);
    }

    public static implicit operator double[,](Matrix m)
    {
        return (double[,])m._values.Clone();
    }

    private static void ValidateSameDimensions(Matrix a, Matrix b)
    {
        if (a.M != b.M || a.N != b.N)
        {
            throw new ArgumentException("Matrix have different dimensions");
        }
    }

    private void ValidateRowIndex(int row)
    {
        if (row < 0 || row >= M)
        {
            throw new ArgumentException("Row index is out of bounds");
        }
    }

    private void ValidateColumnIndex(int column)
    {
        if (column < 0 || column >= N)
        {
            throw new ArgumentException("Column index is out of bounds");
        }
    }

    private static Matrix Operation(Matrix a, Matrix b, 
        Func<double, double, double> operation)
    {
        var result = new double[a.M, a.N];

        for (int i = 0; i < a.M; i++)
        {
            for (int j = 0; j < a.N; j++)
            {
                result[i, j] = operation(a[i, j], b[i, j]);
            }
        }

        return result;
    }

    public static Matrix GetIdentityMatrix(int order)
    {
        var identity = new double[order, order];
        for (int i = 0; i < order; i++)
        {
            for (int j = 0; j < order; j++)
            {
                if (i == j)
                {
                    identity[i, j] = 1;
                }
            }
        }

        return identity;
    }

    public override string ToString()
    {
        var sb = new StringBuilder();
        for (int i = 0; i < M; i++)
        {
            for (int j = 0; j < N; j++)
            {
                sb.Append(string.Format("{0,-8:0.###} ", this[i, j]));
            }
            sb.AppendLine();
        }

        return sb.ToString();
    }

    public override bool Equals(object obj)
    {
        if (ReferenceEquals(null, obj)) return false;
        if (ReferenceEquals(this, obj)) return true;
        if (obj.GetType() != this.GetType()) return false;

        return Equals((Matrix) obj);
    }

    protected bool Equals(Matrix other)
    {
        return _values.Cast<double>()
            .SequenceEqual(other._values.Cast<double>())
            && M == other.M && N == other.N;
    }

    public override int GetHashCode()
    {
        unchecked
        {
            int hashCode = (_values != null ? _values.GetHashCode() : 0);
            hashCode = (hashCode * 397) ^ M;
            hashCode = (hashCode * 397) ^ N;
            return hashCode;
        }
    }
}

Answer 1

It would be better to do the faster checks first in this method:

protected bool Equals(Matrix other)
{
    return _values.Cast<double>()
        .SequenceEqual(other._values.Cast<double>())
        && M == other.M && N == other.N;
}

So we have this:

protected bool Equals(Matrix other)
{
    return M == other.M &&
           N == other.N &&
          _values.Cast<double>().SequenceEqual(other._values.Cast<double>());
}

Cast<double>.SequenceEqual is nice and short, but there might be performance issues (you would need to test this). I think the more verbose option would be better in terms of clarity (and possibly performance):

if (M != other.M || N != other.N)
{
    return false;
}

for (int i = 0; i < M; i++)
{
    for (int j = 0; j < N; j++)
    {
        if (_values[i, j] != other._values[i, j])
        {
            return false;
        }
    }
}

return true;

There's no reason for Equals(Matrix) to be protected; in fact, if we make it public we can let Matrix implement IEquatable<Matrix> (not forgetting to add a null check).

---

On comparing doubles, I'll quote from MSDN:

The Equals method should be used with caution, because two apparently equivalent values can be unequal due to the differing precision of the two values. The following example reports that the Double value .333333 and the Double value returned by dividing 1 by 3 are unequal.

...

Rather than comparing for equality, one technique involves defining an acceptable relative margin of difference between two values (such as .001% of one of the values). If the absolute value of the difference between the two values is less than or equal to that margin, the difference is likely to be due to differences in precision and, therefore, the values are likely to be equal.

---

public static Matrix GetIdentityMatrix(int order)
{
    var identity = new double[order, order];
    for (int i = 0; i < order; i++)
    {
        for (int j = 0; j < order; j++)
        {
            if (i == j)
            {
                identity[i, j] = 1;
            }
        }
    }

    return identity;
}

We can write this in a more efficient way:

public static Matrix GetIdentityMatrix(int order)
{
    var identity = new double[order, order];
    for (int i = 0; i < order; i++)
    {
        identity[i, i] = 1;
    }

    return identity;
}

We should also be checking that order is positive.