Eigen value calculations inconsistent between literals/non-literals and datatypes #1126
Description
WHEN Calculating eigen values from a DenseMatrix or Matrix<T> populated with the same numeric values
GIVEN the numeric datatypes present in DenseMatrix or Matrix<T> (i.e. double/float vs only double vs only float)
GIVEN the calculation contains a non-literal number vs a literal estimate of the number
THEN the output of the eigen value calculation can be different
SHOULD output should be equal assuming all numeric values are equivalent
In human speak:
The output when calculating eigen values from a DenseMatrix or Matrix<T> can change given factors that mathematically shouldn't impact the output. These are:
Matrix<double>vsMatrix<float>vs DenseMatrix (only populated with one data type) vs DenseMatrix (populated with both doubles and floats)- whether the calculation contains a literal or non-literal number
The numeric values the data represent may be the same but the output can differ.
This behaviour is not consistently observable and requires the matrix to undergo transformation processes to trigger. An example is provided below of how this is triggered.
See the screenshot and code snippet below for the example.
Matrix A: DenseMatrix(doubles) + double non-literal
Matrix B: DenseMatrix(doubles) + float non-literal
Matrix C: Matrix<double> + double non-literal
Matrix D: DenseMatrix(doubles) + double literal estimate
Matrix E: DenseMatrix(doubles) + float literal estimate
Matrix F:Matrix<double>+ double literal estimate
Matrix G: Matrix<float> + float non-literal
Matrix H: Matrix<float> + float literal estimate
The offending line(s) that triggers the behavioural changes are when we assign element [1,1] in each matrix which is the only difference in assigning values.
double literalD = 1.96078; //~ 1/0.51
float literalF = 1.96078f; //~ 1/0.51
a[1, 1] = 1 / 0.51; //double non-literal
b[1, 1] = 1 / 0.51f; //float non-literal
c[1, 1] = 1 / 0.51; //double non-literal
d[1, 1] = literalD; //double
e[1, 1] = literalF; //float
f[1, 1] = literalD; //double
g[1, 1] = 1 / 0.51f; //float non-literal
h[1, 1] = literalF; //float
Several pre-processing calculations were done in the scenario that initially triggered this and are required to trigger it here too. When observing in a debugger the final state of these matrices before doing the eigen value calculations the numeric values and indexes were equal.
Out of the above we end up with 2 distinct groups:
A,C,E,H all behave the same.
B,C,F,G all behave the same.
The magnitude of the output values are the same/within the expected deviation given significant figures. However the index of the values within the output matrices is different and sometimes the polarity of the elements flipped.
[TestMethod]
public void FooTest()
{
//define several arrangements of matrices to do operations on.
DenseMatrix a = new DenseMatrix(2, 2); //All doubles - contains non-literal numeral
DenseMatrix b = new DenseMatrix(2, 2); //Mixed doubles/floats - contains non-literal numeral
Matrix<double> c = Matrix<double>.Build.Dense(2, 2); //All doubles - contains non-literal numeral
DenseMatrix d = new DenseMatrix(2, 2); //All doubles - contains literal estimate of value
DenseMatrix e = new DenseMatrix(2, 2); //Mixed doubles/floats - contains literal estimate of value
Matrix<double> f = Matrix<double>.Build.Dense(2, 2); //All doubles - contains literal estimate of value
Matrix<float> g = Matrix<float>.Build.Dense(2, 2); //All floats - contains non-literal numeral
Matrix<float> h = Matrix<float>.Build.Dense(2, 2); //All floats - contains literal estimate of value
//populate with the same values, the only exception is numeric typing on one field.
//Arbitrary equivalent values
//---------
a[0, 0] = 0.51;
b[0, 0] = 0.51;
c[0, 0] = 0.51;
d[0, 0] = 0.51;
e[0, 0] = 0.51;
f[0, 0] = 0.51;
g[0, 0] = 0.51f;
h[0, 0] = 0.51f;
a[0, 1] = -1.52;
b[0, 1] = -1.52;
c[0, 1] = -1.52;
d[0, 1] = -1.52;
e[0, 1] = -1.52;
f[0, 1] = -1.52;
g[0, 1] = -1.52f;
h[0, 1] = -1.52f;
a[1, 0] = -1.52;
b[1, 0] = -1.52;
c[1, 0] = -1.52;
d[1, 0] = -1.52;
e[1, 0] = -1.52;
f[1, 0] = -1.52;
g[1, 0] = -1.52f;
h[1, 0] = -1.52f;
//----------
//In the final field put the same numeral in different forms
//define literals
double literalD = 1.96078; //~ 1/0.51
float literalF = 1.96078f; //~ 1/0.51
a[1, 1] = 1 / 0.51; //double non-literal
b[1, 1] = 1 / 0.51f; //float non-literal
c[1, 1] = 1 / 0.51; //double non-literal
d[1, 1] = literalD; //double
e[1, 1] = literalF; //float
f[1, 1] = literalD; //double
g[1, 1] = 1 / 0.51f; //float non-literal
h[1, 1] = literalF; //float
//Do the same pre-processing calculations
var aT = a.TransposeAndMultiply(a);
var bT = b.TransposeAndMultiply(b);
var cT = c.TransposeAndMultiply(c);
var dT = d.TransposeAndMultiply(d);
var eT = e.TransposeAndMultiply(e);
var fT = f.TransposeAndMultiply(f);
var gT = g.TransposeAndMultiply(g);
var hT = h.TransposeAndMultiply(h);
var aTI = aT.Inverse();
var bTI = bT.Inverse();
var cTI = cT.Inverse();
var dTI = dT.Inverse();
var eTI = eT.Inverse();
var fTI = fT.Inverse();
var gTI = gT.Inverse();
var hTI = hT.Inverse();
//Calculate eigen values.
var aEvd = aTI.Evd();
var bEvd = bTI.Evd();
var cEvd = cTI.Evd();
var dEvd = dTI.Evd();
var eEvd = eTI.Evd();
var fEvd = fTI.Evd();
var gEvd = gTI.Evd();
var hEvd = hTI.Evd();
//Compare the output
Console.WriteLine($"eigen values: a[0] = {aEvd.EigenValues[0]}, b[0] = {bEvd.EigenValues[0]}, c[0] = {cEvd.EigenValues[0]}, d[0] = {dEvd.EigenValues[0]}, " +
$"e[0] = {eEvd.EigenValues[0]}, f[0] = {fEvd.EigenValues[0]}, g[0] = {gEvd.EigenValues[0]}, h[0] = {hEvd.EigenValues[0]}");
Console.WriteLine($"eigen values: a[1] = {aEvd.EigenValues[1]}, b[1] = {bEvd.EigenValues[1]}, c[1] = {cEvd.EigenValues[1]}, d[1] = {dEvd.EigenValues[1]}, " +
$"e[1] = {eEvd.EigenValues[1]}, f[1] = {fEvd.EigenValues[1]}, g[1] = {gEvd.EigenValues[1]}, h[1] = {hEvd.EigenValues[1]}");
Console.WriteLine($"eigen vector: a[0,0] = {aEvd.EigenVectors[0, 0]}, b[0,0] = {bEvd.EigenVectors[0, 0]}, c[0,0] = {cEvd.EigenVectors[0, 0]}, d[0,0] = {dEvd.EigenVectors[0, 0]}, " +
$"e[0,0] = {eEvd.EigenVectors[0, 0]}, f[0,0] = {fEvd.EigenVectors[0, 0]}, g[0,0] = {gEvd.EigenVectors[0, 0]}, h[0,0] = {hEvd.EigenVectors[0, 0]}");
Console.WriteLine($"eigen vector: a[0,1] = {aEvd.EigenVectors[0, 1]}, b[0,1] = {bEvd.EigenVectors[0, 1]}, c[0,1] = {cEvd.EigenVectors[0, 1]}, d[0,1] = {dEvd.EigenVectors[0, 1]}, " +
$"e[0,1] = {eEvd.EigenVectors[0, 1]}, f[0,1] = {fEvd.EigenVectors[0, 1]}, g[0,1] = {gEvd.EigenVectors[0, 1]}, h[0,1] = {hEvd.EigenVectors[0, 1]}");
Console.WriteLine($"eigen vector: a[1,0] = {aEvd.EigenVectors[1, 0]}, b[1,0] = {bEvd.EigenVectors[1, 0]}, c[1,0] = {cEvd.EigenVectors[1, 0]}, d[1,0] = {dEvd.EigenVectors[1, 0]}, " +
$"e[1,0] = {eEvd.EigenVectors[1, 0]}, f[1,0] = {fEvd.EigenVectors[1, 0]}, g[1,0] = {gEvd.EigenVectors[1, 0]}, h[1,0] = {hEvd.EigenVectors[1, 0]}");
Console.WriteLine($"eigen vector: a[1,1] = {aEvd.EigenVectors[1, 1]}, b[1,1] = {bEvd.EigenVectors[1, 1]}, c[1,1] = {cEvd.EigenVectors[1, 1]}, d[1,1] = {dEvd.EigenVectors[1, 1]}, " +
$"e[1,1] = {eEvd.EigenVectors[1, 1]}, f[1,1] = {fEvd.EigenVectors[1, 1]}, g[1,1] = {gEvd.EigenVectors[1, 1]}, h[1,1] = {hEvd.EigenVectors[1, 1]}");
}
output from the above formatted and colour coded to hopefully make the discrepancies more visually obvious
