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

[Martin Maechler]

>>>>> 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.