forum.alglib.net

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;