[Rd] Concerns with SVD -- and the Matrix Exponential
Durga Prasad G me14d059
me14d059 @end|ng |rom @m@||@||tm@@c@|n
Wed Aug 16 10:06:10 CEST 2023
Dear Martin, I am getting different responses from different officials of
R-Software, but those statements are contradicting with the statements
discussed in your email. Kindly go through the previous files and emails,
and respond. I personally think, together we can fix the issue which is
observed in SVD.
Thanks and regards
Durga Prasad
On Tue, Aug 1, 2023 at 4:51 PM Lakshman, Aidan H <AHL27 using pitt.edu> wrote:
> Hi Durga,
>
> There’s an error in your calculations here. You mention that for the SVD
> of a symmetric matrix, we must have U=V, but this is not a correct
> statement. The unitary matrices are only equivalent if the matrix A is
> positive semidefinite.
>
> In your example, you provide the matrix {{1,4},{4,1}}, which has
> eigenvalues 5 and -3. This is not positive semidefinite, and thus there's
> no requirement that the unitary matrices be equivalent.
>
> If you verify your example with something like wolfram alpha, you’ll find
> that R’s solution is correct.
>
> -Aidan
>
> -----------------------
>
> Aidan Lakshman (he/him) <https://www.ahl27.com/>
>
> Doctoral Fellow, Wright Lab <https://www.wrightlabscience.com/>
>
> University of Pittsburgh School of Medicine
>
> Department of Biomedical Informatics
>
> ahl27 using pitt.edu
>
> (724) 612-9940
>
>
>
> ------------------------------
> *From:* R-devel <r-devel-bounces using r-project.org> on behalf of Durga Prasad
> G me14d059 <me14d059 using smail.iitm.ac.in>
> *Sent:* Tuesday, August 1, 2023 4:18:20 AM
> *To:* Martin Maechler <maechler using stat.math.ethz.ch>; r-devel using r-project.org
> <r-devel using r-project.org>; profjcnash using gmail.com <profjcnash using gmail.com>
> *Subject:* Re: [Rd] Concerns with SVD -- and the Matrix Exponential
>
> 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://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdoi.org%2F10.1137%2F1020098&data=05%7C01%7Cahl27%40pitt.edu%7C8575b77db32345ca544b08db927ceae0%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638264837816871329%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=Y4mlFL%2FggLKd7FoIoY62esiFGUwukRG0YmELsJj7nd0%3D&reserved=0
> <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://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdoi.org%2F10.1137%2FS00361445024180&data=05%7C01%7Cahl27%40pitt.edu%7C8575b77db32345ca544b08db927ceae0%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638264837817183809%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=5%2FlssUGC6q7SUy0PY7gZCqWi0%2BXbNwZD0FaAgIcOWdY%3D&reserved=0
> <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://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fcran.r-project.org%2Fpackage%3Dexpm&data=05%7C01%7Cahl27%40pitt.edu%7C8575b77db32345ca544b08db927ceae0%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638264837817183809%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=LlAzmeRZd5tNqAgTuYTTSsSakWEj85%2B%2F%2FP%2FM0DnZNLk%3D&reserved=0
> <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://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.rdocumentation.org%2Fpackages%2Fexpm%2Fversions%2F0.999-7%2Ftopics%2Fexpm&data=05%7C01%7Cahl27%40pitt.edu%7C8575b77db32345ca544b08db927ceae0%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638264837817340048%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=ZrT%2FYvklccqLCORBMf6nGop5o3n7O2thkknG1UAS2Gc%3D&reserved=0
> <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.
> >
>
[[alternative HTML version deleted]]
More information about the R-devel
mailing list