R devel - Aug 2023 - 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.

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