[Rd] Concerns with SVD -- and the Matrix Exponential

Wed Aug 16 18:23:40 CEST 2023

You wrote:
         Using singular value decomposition, any second-order tensor is
given as
          A = UΣVt
          where U and V are the orthogonal tensors, and Σ is the diagonal
matrix (Eigenvalue matrix).

          For a symmetric matrix, the orthogonal tensors are the same,
i.e., U=V.

Can you state your definition of the SVD and prove (or outline a proof of)
that last statement?

-Bill

On Wed, Aug 16, 2023 at 3:47 AM Durga Prasad G me14d059 <
me14d059 using smail.iitm.ac.in> wrote:

> 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]]
>
> ______________________________________________
> R-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

	[[alternative HTML version deleted]]