forum.alglib.net http://forum.alglib.net/ |
|
SVD Accuracy http://forum.alglib.net/viewtopic.php?f=2&t=140 |
Page 1 of 1 |
Author: | LeesNav [ Mon Dec 13, 2010 11:41 pm ] |
Post subject: | SVD Accuracy |
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; |
Author: | Sergey.Bochkanov [ Wed Dec 15, 2010 9:30 am ] |
Post subject: | Re: SVD Accuracy |
There are several reasons behind it: 1. ALGLIB SVD is tuned towards better speed, although you can change it and add about three digits of accuracy. rmatrixsvd() calls rmatrixbdsvd() internally. 5th parameter of rmatrixbdsvd() is false by default; you can set it to true to enforce better precision 2. algorithmic reasons - ALGLIB includes EVD/SVD solvers from LAPACK, but only QR-based ones which are slower and less precise than their modern coutnerparts. It is planned to add divide-and-conquer based EVD/SVD somewhere in 2011 (in the second half, I think) |
Author: | LeesNav [ Wed Dec 15, 2010 10:57 pm ] |
Post subject: | Re: SVD Accuracy |
Thanks for the reply. I definitely appreciate your tuning tips, and I am looking forward to testing the results of divide and conquer based EVD/SVD when they are available. I did a preliminary front to back analysis of my precision requirements, and it looks like the current alglib implementation will get the job done - and at about five times the speed of the dotnumerics solution. So, nice job, and thanks again. |
Page 1 of 1 | All times are UTC |
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group http://www.phpbb.com/ |