First of all, thankyou for making such an excellent library.
I've run into a performance problem with using the C# SVD to implement a Moore-Penrose pseudoinverse. The algorithm works fine when used on smaller matrices, but on larger matrices (12x12) the result is not good. By "not good" I mean that it does not agree well with test vectors created using MATLAB. When I implement the Moore-Penrose pseudoinverse using the DotNumerics implementation of the SVD, the output agrees with the MATLAB test vectors. However, alglib is much faster than DotNumerics. If the accuracy of the alglib svd was improved, that would be awesome.
I posted a large chunk of code to help demonstrate what I am referring to. I hope it makes things clear:
It contains:
the input test vector
the Moore-Penrose pseudoinverse implementation
the expected output (generated from MATLAB.)
Code:
double[,] input_matrix = new double[12,12];
input_matrix[0, 0] = 0.27958153633994054;
input_matrix[0, 1] = -0.0026248302281403971;
input_matrix[0, 2] = 0.091102023161413914;
input_matrix[0, 3] = 0.020502256732689038;
input_matrix[0, 4] = 0.0;
input_matrix[0, 5] = 0.0;
input_matrix[0, 6] = 0.0;
input_matrix[0, 7] = 0.0;
input_matrix[0, 8] = -152.4238315838742;
input_matrix[0, 9] = 38.037963410472784;
input_matrix[0, 10] = -640.72479822146863;
input_matrix[0, 11] = -117.02403520736323;
input_matrix[1, 0] = -0.0026248302281403971;
input_matrix[1, 1] = 0.24055788933609173;
input_matrix[1, 2] = -0.56710135350560109;
input_matrix[1, 3] = -0.089817031070816483;
input_matrix[1, 4] = 0.0;
input_matrix[1, 5] = 0.0;
input_matrix[1, 6] = 0.0;
input_matrix[1, 7] = 0.0;
input_matrix[1, 8] = 38.0379634104728;
input_matrix[1, 9] = -114.24283232651035;
input_matrix[1, 10] = 275.03889446852617;
input_matrix[1, 11] = 43.322683830078269;
input_matrix[2, 0] = 0.091102023161413914;
input_matrix[2, 1] = -0.56710135350560109;
input_matrix[2, 2] = 13.999999999999996;
input_matrix[2, 3] = 2.4588664235958366;
input_matrix[2, 4] = 0.0;
input_matrix[2, 5] = 0.0;
input_matrix[2, 6] = 0.0;
input_matrix[2, 7] = 0.0;
input_matrix[2, 8] = -640.72479822146852;
input_matrix[2, 9] = 275.03889446852611;
input_matrix[2, 10] = -6846.0773607294814;
input_matrix[2, 11] = -1212.0190360038757;
input_matrix[3, 0] = 0.020502256732689038;
input_matrix[3, 1] = -0.089817031070816483;
input_matrix[3, 2] = 2.4588664235958366;
input_matrix[3, 3] = 0.43412934612588078;
input_matrix[3, 4] = 0.0;
input_matrix[3, 5] = 0.0;
input_matrix[3, 6] = 0.0;
input_matrix[3, 7] = 0.0;
input_matrix[3, 8] = -117.02403520736323;
input_matrix[3, 9] = 43.322683830078248;
input_matrix[3, 10] = -1212.019036003876;
input_matrix[3, 11] = -215.98768637506004;
input_matrix[4, 0] = 0.0;
input_matrix[4, 1] = 0.0;
input_matrix[4, 2] = 0.0;
input_matrix[4, 3] = 0.0;
input_matrix[4, 4] = 0.27958153633994054;
input_matrix[4, 5] = -0.0026248302281403971;
input_matrix[4, 6] = 0.091102023161413914;
input_matrix[4, 7] = 0.020502256732689038;
input_matrix[4, 8] = -119.9223539336146;
input_matrix[4, 9] = 1.5168332144743655;
input_matrix[4, 10] = -45.445137604688043;
input_matrix[4, 11] = -10.837929007543186;
input_matrix[5, 0] = 0.0;
input_matrix[5, 1] = 0.0;
input_matrix[5, 2] = 0.0;
input_matrix[5, 3] = 0.0;
input_matrix[5, 4] = -0.0026248302281403971;
input_matrix[5, 5] = 0.24055788933609173;
input_matrix[5, 6] = -0.56710135350560109;
input_matrix[5, 7] = -0.089817031070816483;
input_matrix[5, 8] = 1.5168332144743637;
input_matrix[5, 9] = -88.8860048198343;
input_matrix[5, 10] = -195.41469566219828;
input_matrix[5, 11] = -34.544267185097922;
input_matrix[6, 0] = 0.0;
input_matrix[6, 1] = 0.0;
input_matrix[6, 2] = 0.0;
input_matrix[6, 3] = 0.0;
input_matrix[6, 4] = 0.091102023161413914;
input_matrix[6, 5] = -0.56710135350560109;
input_matrix[6, 6] = 13.999999999999996;
input_matrix[6, 7] = 2.4588664235958366;
input_matrix[6, 8] = -45.445137604688043;
input_matrix[6, 9] = -195.41469566219834;
input_matrix[6, 10] = -6634.4511546204521;
input_matrix[6, 11] = -1186.1038617201416;
input_matrix[7, 0] = 0.0;
input_matrix[7, 1] = 0.0;
input_matrix[7, 2] = 0.0;
input_matrix[7, 3] = 0.0;
input_matrix[7, 4] = 0.020502256732689038;
input_matrix[7, 5] = -0.089817031070816483;
input_matrix[7, 6] = 2.4588664235958366;
input_matrix[7, 7] = 0.43412934612588078;
input_matrix[7, 8] = -10.837929007543186;
input_matrix[7, 9] = -34.544267185097908;
input_matrix[7, 10] = -1186.1038617201416;
input_matrix[7, 11] = -213.11562747737682;
input_matrix[8, 0] = -152.4238315838742;
input_matrix[8, 1] = 38.0379634104728;
input_matrix[8, 2] = -640.72479822146852;
input_matrix[8, 3] = -117.02403520736323;
input_matrix[8, 4] = -119.9223539336146;
input_matrix[8, 5] = 1.5168332144743637;
input_matrix[8, 6] = -45.445137604688043;
input_matrix[8, 7] = -10.837929007543186;
input_matrix[8, 8] = 189833.58531503039;
input_matrix[8, 9] = -37489.998041992774;
input_matrix[8, 10] = 652352.43902029027;
input_matrix[8, 11] = 120626.89004661031;
input_matrix[9, 0] = 38.037963410472784;
input_matrix[9, 1] = -114.24283232651035;
input_matrix[9, 2] = 275.03889446852611;
input_matrix[9, 3] = 43.322683830078248;
input_matrix[9, 4] = 1.5168332144743655;
input_matrix[9, 5] = -88.8860048198343;
input_matrix[9, 6] = -195.41469566219834;
input_matrix[9, 7] = -34.544267185097908;
input_matrix[9, 8] = -37489.998041992774;
input_matrix[9, 9] = 125875.22355671997;
input_matrix[9, 10] = 87148.863552965471;
input_matrix[9, 11] = 17359.754529162197;
input_matrix[10, 0] = -640.72479822146863;
input_matrix[10, 1] = 275.03889446852617;
input_matrix[10, 2] = -6846.0773607294814;
input_matrix[10, 3] = -1212.019036003876;
input_matrix[10, 4] = -45.445137604688043;
input_matrix[10, 5] = -195.41469566219828;
input_matrix[10, 6] = -6634.4511546204521;
input_matrix[10, 7] = -1186.1038617201416;
input_matrix[10, 8] = 652352.43902029027;
input_matrix[10, 9] = 87148.863552965471;
input_matrix[10, 10] = 8755890.8047270589;
input_matrix[10, 11] = 1564165.090420241;
input_matrix[11, 0] = -117.02403520736323;
input_matrix[11, 1] = 43.322683830078269;
input_matrix[11, 2] = -1212.0190360038757;
input_matrix[11, 3] = -215.98768637506004;
input_matrix[11, 4] = -10.837929007543186;
input_matrix[11, 5] = -34.544267185097922;
input_matrix[11, 6] = -1186.1038617201416;
input_matrix[11, 7] = -213.11562747737682;
input_matrix[11, 8] = 120626.89004661031;
input_matrix[11, 9] = 17359.754529162197;
input_matrix[11, 10] = 1564165.090420241;
input_matrix[11, 11] = 281146.31885908241;
///////////////////////////////////////////////////////////////
//
// Moore-Penrose pseudoinverse implementation
//
///////////////////////////////////////////////////////////////
double[] singular_values = new double[input_matrix.GetLength(1)];
double[,] u = new double[input_matrix.GetLength(0), input_matrix.GetLength(0)];
double[,] v_transpose = new double[input_matrix.GetLength(1), input_matrix.GetLength(1)];
double[,] singular_values_reciprocal = new double[input_matrix.GetLength(0), input_matrix.GetLength(1)];
double[,] temp_matrix = new double[input_matrix.GetLength(1), input_matrix.GetLength(0)];
double[,] output_matrix = new double[input_matrix.GetLength(1), input_matrix.GetLength(0)];
alglib.svd.rmatrixsvd(
input_matrix, // input
input_matrix.GetLength(0), // number of rows
input_matrix.GetLength(1), // number of columns
2, // compute the whole of U
2, // compute the whole of V_transpose
2, // use maximum memory for more performance
ref singular_values, // output singular values, m length vector
ref u, // output u, m by m matrix
ref v_transpose); // output v_transpose, n by n matrix
// Compute the smallest value that will be considered in computing the reciprocal of r1_singular_values
// smaller non zero values that should be zero may occur due to numerical precision issues.
// for example, the reciprocal of 1e-17 would be huge, and 1e-17 is at the limit of double precision.
double largest_singular_value = singular_values[0];
double max_dimension = 0;
if (input_matrix.GetLength(0) > input_matrix.GetLength(1))
{
max_dimension = input_matrix.GetLength(0);
}
else
{
max_dimension = input_matrix.GetLength(1);
}
double machine_precision_limit = Math.Pow(1.11 * 10, -16);
double minimum_nonzero_value = largest_singular_value * max_dimension * machine_precision_limit;
// Compute the matrix composed of reciprocals of the r1_singular_values vector
for (uint row_index = 0; row_index < input_matrix.GetLength(0); row_index++)
{
for (uint column_index = 0; column_index < input_matrix.GetLength(1); column_index++)
{
if (column_index == row_index)
{
if (Math.Abs(singular_values[column_index]) > minimum_nonzero_value)
{
singular_values_reciprocal[row_index, column_index] = 1.0 / singular_values[column_index];
}
else
{
singular_values_reciprocal[row_index, column_index] = 0.0;
}
}
else
{
singular_values_reciprocal[row_index, column_index] = 0.0;
}
}
}
// Multiply transpose(r2_v_transpose) and transpose(r2_singular_values_reciprocal)
ablas.rmatrixgemm(
input_matrix.GetLength(1),
input_matrix.GetLength(0),
input_matrix.GetLength(1),
1.0,
ref v_transpose, // n x n matrix
0,
0,
1,
ref singular_values_reciprocal, // m x n matrix
0,
0,
1, // take the transpose of singular_values_reciprocal
0,
ref temp_matrix, // n x m matrix
0,
0);
// Multiply temporary_matrix_one and transpose(u)
ablas.rmatrixgemm(
input_matrix.GetLength(1),
input_matrix.GetLength(0),
input_matrix.GetLength(0),
1.0,
ref temp_matrix, // n x m matrix
0,
0,
0,
ref u, // m x m matrix
0,
0,
1, // take the transpose of u
0,
ref output_matrix,
0,
0);
// output_matrix contains the Moore-Penrose pseudoinverse
// Now I will list the expected result, as produced by MATLAB. Note that the pseudoinverse
// implemented using the dotnumerics svd matches these results to 10^-7.
double[,]expected = new double[12, 12];
expected[0, 0] = 204.53347710446158;
expected[0, 1] = -146.97687695504621;
expected[0, 2] = -183.66036583299717;
expected[0, 3] = 1229.5992242275388;
expected[0, 4] = 3.2817707556157427;
expected[0, 5] = -79.144734796839487;
expected[0, 6] = -323.00722339668243;
expected[0, 7] = 1986.7527225703329;
expected[0, 8] = 0.02525371577531214;
expected[0, 9] = -0.30033110018715903;
expected[0, 10] = -0.3090509536197944;
expected[0, 11] = 2.1214747491021919;
expected[1, 0] = -146.97685293235446;
expected[1, 1] = 115.74829467162877;
expected[1, 2] = 129.78549441887534;
expected[1, 3] = -871.53988515374442;
expected[1, 4] = -1.0413588895094523;
expected[1, 5] = 58.050161824389157;
expected[1, 6] = 229.48872429565998;
expected[1, 7] = -1413.6240541088941;
expected[1, 8] = -0.01507424860467569;
expected[1, 9] = 0.22718537482569126;
expected[1, 10] = 0.21492992143413389;
expected[1, 11] = -1.4886876346382281;
expected[2, 0] = -183.66032537796951;
expected[2, 1] = 129.78548914214392;
expected[2, 2] = 207.09217315822269;
expected[2, 3] = -1349.2713496931983;
expected[2, 4] = -2.8315159996070474;
expected[2, 5] = 58.8795718839225;
expected[2, 6] = 300.68580282194722;
expected[2, 7] = -1856.8088731500998;
expected[2, 8] = -0.020270805856085389;
expected[2, 9] = 0.25395677425834046;
expected[2, 10] = 0.31628682191877988;
expected[2, 11] = -2.138733426397025;
expected[3, 0] = 1229.5990051681201;
expected[3, 1] = -871.53988550348788;
expected[3, 2] = -1349.2713984041386;
expected[3, 3] = 8828.362807710193;
expected[3, 4] = 13.5612279594478;
expected[3, 5] = -402.85940008897825;
expected[3, 6] = -2011.6441498501044;
expected[3, 7] = 12423.590415802317;
expected[3, 8] = 0.1255997716194508;
expected[3, 9] = -1.7212956783708504;
expected[3, 10] = -2.0894119613271074;
expected[3, 11] = 14.170250613252382;
expected[4, 0] = 3.2817708588782462;
expected[4, 1] = -1.0413588918952978;
expected[4, 2] = -2.8315159788440254;
expected[4, 3] = 13.561227971045065;
expected[4, 4] = 10.810903199043299;
expected[4, 5] = -3.1428285510185261;
expected[4, 6] = 1.8502794642060914;
expected[4, 7] = -17.389446526260954;
expected[4, 8] = 0.016730054569402741;
expected[4, 9] = 0.00219956977460149;
expected[4, 10] = 0.0013544520122326;
expected[4, 11] = -0.0204563787374225;
expected[5, 0] = -79.144689274575853;
expected[5, 1] = 58.050137902193605;
expected[5, 2] = 58.87954158920008;
expected[5, 3] = -402.85918123488813;
expected[5, 4] = -3.1428286507134842;
expected[5, 5] = 35.49616006160894;
expected[5, 6] = 116.77833357507861;
expected[5, 7] = -715.45845919694989;
expected[5, 8] = -0.01420756877953187;
expected[5, 9] = 0.120151136166711;
expected[5, 10] = 0.10447235913425988;
expected[5, 11] = -0.72553762106404551;
expected[6, 0] = -323.00715910966062;
expected[6, 1] = 229.48871809910085;
expected[6, 2] = 300.68580650387946;
expected[6, 3] = -2011.6440923810044;
expected[6, 4] = 1.8502794410497005;
expected[6, 5] = 116.77838585980382;
expected[6, 6] = 532.33928331511;
expected[6, 7] = -3278.8221855436268;
expected[6, 8] = -0.025034079295335129;
expected[6, 9] = 0.472637610568542;
expected[6, 10] = 0.51273521584087656;
expected[6, 11] = -3.5152059260915722;
expected[7, 0] = 1986.7523806556278;
expected[7, 1] = -1413.6240546565009;
expected[7, 2] = -1856.8089491616593;
expected[7, 3] = 12423.590415705678;
expected[7, 4] = -17.389446546697911;
expected[7, 5] = -715.45880077254446;
expected[7, 6] = -3278.8222752145371;
expected[7, 7] = 20207.800938335302;
expected[7, 8] = 0.13980711139461249;
expected[7, 9] = -2.9125793310887369;
expected[7, 10] = -3.1635421814857008;
expected[7, 11] = 21.701386336428808;
expected[8, 0] = 0.025253713089920221;
expected[8, 1] = -0.015074248614759429;
expected[8, 2] = -0.020270806457849282;
expected[8, 3] = 0.12559977164386851;
expected[8, 4] = 0.01673005456919863;
expected[8, 5] = -0.014207571469796241;
expected[8, 6] = -0.025034080008072751;
expected[8, 7] = 0.13980711143868205;
expected[8, 8] = 0.0000409790721521;
expected[8, 9] = -0.00002097066783035;
expected[8, 10] = -0.00002027022042745;
expected[8, 11] = 0.00011768668472994;
expected[9, 0] = -0.3003310513068751;
expected[9, 1] = 0.22718537485757329;
expected[9, 2] = 0.25395678504531904;
expected[9, 3] = -1.7212956779285342;
expected[9, 4] = 0.00219956977884596;
expected[9, 5] = 0.12015118490353421;
expected[9, 6] = 0.47263762325872755;
expected[9, 7] = -2.9125793304059262;
expected[9, 8] = -0.00002097066781507;
expected[9, 9] = 0.00047214678544382;
expected[9, 10] = 0.00044101842814172;
expected[9, 11] = -0.00306033986798547;
expected[10, 0] = -0.30905088327916186;
expected[10, 1] = 0.21492991031975972;
expected[10, 2] = 0.31628681864907038;
expected[10, 3] = -2.0894118589255517;
expected[10, 4] = 0.00135445196763594;
expected[10, 5] = 0.10447240813192527;
expected[10, 6] = 0.51273520424409091;
expected[10, 7] = -3.1635420216909496;
expected[10, 8] = -0.00002027021916462;
expected[10, 9] = 0.00044101840543863;
expected[10, 10] = 0.00053629123965554;
expected[10, 11] = -0.00362763995095133;
expected[11, 0] = 2.1214743759367587;
expected[11, 1] = -1.4886876353134946;
expected[11, 2] = -2.1387335094946125;
expected[11, 3] = 14.170250613892215;
expected[11, 4] = -0.02045637875668533;
expected[11, 5] = -0.725537994024669;
expected[11, 6] = -3.5152060241834344;
expected[11, 7] = 21.701386337617762;
expected[11, 8] = 0.00011768668469893;
expected[11, 9] = -0.00306033986889718;
expected[11, 10] = -0.00362764012559858;
expected[11, 11] = 0.024633259441036559;