Algorithms

public static Matrix CovarianceMatrix(List<double[]> dataset)
        {
            /*
                cov_xx cov_xy ...
                cov_yx cov_yy ...
                ...
             */

            //Calculate mean for this cluster
            // cov_xx = sum[x*x]/n, cov_xy = sum[x*y]/n
            double[] means = new double[dataset[0].Length];
            Matrix cov = new Matrix(dataset[0].Length, dataset[0].Length);
            double sum = 0;

            for (int i = 0; i < dataset[0].Length; i++)
            {
                for (int j = 0; j < dataset.Count; j++)
                {
                    means[i] += dataset[j][i];
                }
                means[i] /= dataset.Count;
            }

            double[,] subresults = new double[dataset[0].Length, dataset.Count];
            for (int j = 0; j < dataset.Count; j++)
            {
                for (int i = 0; i < dataset[0].Length; i++)
                {
                    subresults[i, j] = dataset[j][i] - means[i];
                }
            }
            
            //fill covariance
            for (int i = 0; i < dataset[0].Length; i++)
            {
                for (int j = i; j < dataset[0].Length; j++)
                {
                    double s = 0;
                    for (int x = 0; x < dataset.Count; x++)
                    {
                        s += subresults[i, x] * subresults[j, x];
                    }
                    cov.SetElement(i, j, s / dataset.Count);
                    if (i != j) cov.SetElement(j, i, s / dataset.Count);
                }
            }
            return cov;
        }

public static double Mahalanobis(double[] vector1, double[] vector2, Matrix covariance)
{
    Matrix v1 = new Matrix(vector1, vector1.Length);
    Matrix v2 = new Matrix(vector2, vector2.Length);
    Matrix m = v1.Subtract(v2);

    return (double)(m.Transpose()).Multiply(covariance.Inverse()).Multiply(m).GetElement(0, 0);
}