[Rd] Concerns with SVD -- and the Matrix Exponential
Durga Prasad G me14d059
me14d059 @end|ng |rom @m@||@||tm@@c@|n
Tue Aug 1 10:18:20 CEST 2023
Hi Martin, Thank you for your reply. The response and the links provided by
you helped to learn more. But I am not able to obtain the simple even
powers of a matrix: one simple case is the square of a matrix. The square
of the matrix using direct matrix multiplication operations and svd (A = U
D V') are different. Kindly check the attached file for the complete
explanation. I want to know which technique was used in building the svd in
R-Software. I want to discuss about svd if you schedule a meeting.
Thanks and Regards
Durga Prasad
On Mon, Jul 17, 2023 at 2:13 PM Martin Maechler <maechler using stat.math.ethz.ch>
wrote:
> >>>>> 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.
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Square.pdf
Type: application/pdf
Size: 79039 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20230801/8c106736/attachment.pdf>
More information about the R-devel
mailing list