Better check your definitions of SVD -- there are several forms, but all I am aware of (and I wrote a couple of the codes in the early 1970s for the SVD) have positive singular values. JN On 2023-07-16 02:01, Durga Prasad G me14d059 wrote:> Respected Development Team, > > This is Durga Prasad reaching out to you regarding an issue/concern related > to Singular Value Decomposition SVD in R software package. I am attaching a > detailed attachment with this letter which depicts real issues with SVD in > R. > > To reach the concern the expressions for the exponential of a matrix using > SVD and > projection tensors are obtained from series expansion. However, numerical > inconsistency is observed between the exponential of matrix obtained using > the function(svd()) used in R software. > > However, it is observed that most of the researchers fraternity is engaged > in utilising R software for their research purposes and to the extent of my > understanding such an error in SVD in R software might raise the concern > about authenticity of the simulation results produced and published by > researchers across the globe. > > Further, I am very sure that the R software development team is well versed > with the happening and they have any specific and resilient reasons for > doing so. I would request you kindly, to guide me through the concern. > > Thank you very much. > > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel
Martin Maechler
2023-Jul-17 08:43 UTC
[Rd] Concerns with SVD -- and the Matrix Exponential
>>>>> J C Nash >>>>> on Sun, 16 Jul 2023 13:30:57 -0400 writes:> Better check your definitions of SVD -- there are several > forms, but all I am aware of (and I wrote a couple of the > codes in the early 1970s for the SVD) have positive > singular values. > JN Indeed. More generally, the decomposition A = U D V' (with diagonal D and orthogonal U,V) is not at all unique. There are not only many possible different choices of the sign of the diagonal entries, but also the *ordering* of the singular values is non unique. That's why R and 'Lapack', the world-standard for computer/numerical linear algebra, and others I think, make the decomposition unique by requiring non-negative entries in D and have them *sorted* decreasingly. The latter is what the help page help(svd) always said (and you should have studied that before raising such concerns). ----------------------------------------------------------------- To your second point (in the document), the matrix exponential: It is less known, but still has been known among experts for many years (and I think even among students of a class on numerical linear algebra), that there are quite a few mathematically equivalent ways to compute the matrix exponential, *BUT* that most of these may be numerically disastrous, for several different reasons depending on the case. This has been known for close to 50 years now: Cleve Moler and Charles Van Loan (1978) Nineteen Dubious Ways to Compute the Exponential of a Matrix SIAM Review Vol. 20(4) https://doi.org/10.1137/1020098 Where as that publication had been important and much cited at the time, the same authors (known world experts in the field) wrote a review of that review 25 years later which I think (and hope) is even more widely cited (in R's man/*.Rd syntax) : Cleve Moler and Charles Van Loan (2003) Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. \emph{SIAM Review} \bold{45}, 1, 3--49. \doi{10.1137/S00361445024180} i.e. https://doi.org/10.1137/S00361445024180 It is BTW also cited on the Wikipedia page on the matrix exponential: ==> For this reason, Professor Douglas Bates, the initial creator of R's Matrix package (which comes with R) has added the Matrix exponential very early to the package: ------------------------------------------------------------------------ r461 | bates | 2005-01-29 Add expm function ------------------------------------------------------------------------ Later, I've fixed an "infamous" bug : ------------------------------------------------------------------------ r2127 | maechler | 2008-03-07 fix the infamous expm() bug also in "Matrix" (duh!) ------------------------------------------------------------------------ Then, Vincent Goulet and Christophe Dutang wanted to provide more versions of expm() and we collaborated, also providing expm() for complex matrices and created the CRAN package {expm} --> https://CRAN.R-project.org/package=expm which already provided half a dozen different expm algorithms. But research and algorithms did not stop there. In 2008, Higham and collaborators even improved on the best known algorithms and I had the chance to co-supervise a smart Master's student Michael Stadelmann to implement Higham's algorithm and we even allowed to tweak it {with optional arguments} as that was seen to be beneficial sometimes. See e.g., https://www.rdocumentation.org/packages/expm/versions/0.999-7/topics/expm > On 2023-07-16 02:01, Durga Prasad G me14d059 wrote: >> Respected Development Team, >> >> This is Durga Prasad reaching out to you regarding an >> issue/concern related to Singular Value Decomposition SVD >> in R software package. I am attaching a detailed >> attachment with this letter which depicts real issues >> with SVD in R. >> >> To reach the concern the expressions for the exponential >> of a matrix using SVD and projection tensors are obtained >> from series expansion. However, numerical inconsistency >> is observed between the exponential of matrix obtained >> using the function(svd()) used in R software. >> >> However, it is observed that most of the researchers >> fraternity is engaged in utilising R software for their >> research purposes and to the extent of my understanding >> such an error in SVD in R software might raise the >> concern about authenticity of the simulation results >> produced and published by researchers across the globe. >> >> Further, I am very sure that the R software development >> team is well versed with the happening and they have any >> specific and resilient reasons for doing so. I would >> request you kindly, to guide me through the concern. >> >> Thank you very much.