Thanks for the feedback! It turned out that SVD was fine and I was looking in the wrong place -- in the C counterpart, I was feeding my ellipsoid algorithm with uninitialized junk, and the compiler didn't shout at me.
Sorry for the noise! niedz., 26 wrz 2021 o 15:16 Eli Zaretskii <e...@gnu.org> napisał(a): > > From: Panicz Maciej Godek <godek.mac...@gmail.com> > > Date: Sun, 26 Sep 2021 13:55:20 +0200 > > > > I forgot to mention that I run Guile in an Emacs session running in a WSL > > console on Windows 10. > > The tests of the C code that I've been performing so far were executed in > > an MSYS terminal, but I have just tried running them in WSL console, and > > the radii get ridiculous values. > > > > While the values used to generate the points are the following: > > > > center '[-65.12 -50.54 88.66]: > > radii: '[83.95 47.13 45.56] > > rotation: '[[9.633871e-001 1.363818e-001 -2.308359e-001] > > [-1.734094e-001 9.735887e-001 -1.485064e-001] > > [2.044857e-001 1.830983e-001 9.615927e-001]] > > > > and the values reconstructed in Guile are: > > > > (ellipsoid > > #:center (65.10194623013226 -88.58460232582514 -50.721825868168324) > > #:radii (83.94677717528019 45.56525864722978 47.12037877216948) > > #:rotation ((-0.9633227088329351 0.2307406879308479 > 0.13699669185777974) > > (0.173638969876562 0.14674930301995573 > > 0.9738142277679884) > > (-0.20459439578586436 -0.961885324244193 > > 0.18143251146545056))) > > > > (the signs and the order of radii differ, but it seems that the rotation > > matrix compensates for that) > > > > in the case of MSYS (MinGW 64-bit) I get the following values: > > > > approx_e.center = [90.13, -68.76, -42.46] > > approx_e.radius = [178.83, 97.08, 100.43] > > approx_e.rotation(3x3) = > > -9.634121049254586e-001 1.734020884333972e-001 > > -2.043742444879810e-001 > > -2.314281253565038e-001 -1.535574849495593e-001 > > 9.606566096217422e-001 > > -1.351966674037161e-001 -9.728061546592428e-001 > > -1.880692600613582e-001 > > > > but in the case of the WSL console, I get > > > > approx_e.center = [90.13, -68.76, -42.46] > > approx_e.radius = [1347844355973136335704668472606720.00, > > 731678657375689861259088782950400.00, > 756961648396782369967223865344000.00] > > approx_e.rotation(3x3) = > > -9.634121049254428e-01 1.734020884334219e-01 > > -2.043742444880351e-01 > > -2.314281253566119e-01 -1.535574849499134e-01 > > 9.606566096216593e-01 > > -1.351966674036449e-01 -9.728061546591822e-01 > > -1.880692600617213e-01 > > Forgive me for saying this, but frankly the most probable reason for > this is some bug in converting the Scheme code to C. SVD is a > remarkably stable algorithm in the face of numerical roundoff errors, > that's one of its strong points. So if the issue is with some minor > numerical inaccuracies, SVD results should be insensitive to it, as > long as the condition numbers of the matrices that you decompose are > not preposterously large, like 10^12 or something. (Did you calculate > the condition numbers, btw? what are they?) > > So I would check and recheck your C code, as IMO that's the first > suspect in this story. > > Also, please note that on Intel CPUs double-precision FP calculations > usually keep intermediate results in 10-byte (80-bit) wide "long > double" format, unless your program is built with the GCC switch that > disables that (few programs do, so I don't expect you to do that, > eneither while building Guile, nor when compiling your C > implementation of SVD). But that cannot by itself explain the > problem, because using 80-bit FP values should make the results more > accurate, not less accurate. > > So again, I suspect your C program. Does it work well when you throw > simple problems at it, like systems of linear equations where the > result is known in advance and the condition number of the matrix is > around 1? > > HTH >
