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

2023-08-16 Thread Duncan Murdoch
Dear Durga, I think you have a basic misunderstanding of this mailing 
list.  The responses you have received are from users and volunteer 
developers.  There are no "officials of R-Software".  R is an open 
source project containing contributions from hundreds (maybe thousands) 
of people.


It's only natural that there will be some contradictions in the 
responses from those people.  It's up to you to read the responses and 
find the parts of them that are useful to you.  It's rude of you to ask 
one particular respondent to do that work for you.


Duncan Murdoch

On 16/08/2023 4:06 a.m., Durga Prasad G me14d059 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  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

ah...@pitt.edu

(724) 612-9940



--
*From:* R-devel  on behalf of Durga Prasad
G me14d059 
*Sent:* Tuesday, August 1, 2023 4:18:20 AM
*To:* Martin Maechler ; r-devel@r-project.org
; profjcn...@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 <
maech...@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=05%7C01%7Cahl27%40pitt.edu%7C8575b77db32345ca544b08db927ceae0%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638264837816871329%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=Y4mlFL%2FggLKd7FoIoY62esiFGUwukRG0YmELsJj7nd0%3D=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 syn

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

2023-08-16 Thread Bill Dunlap
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 <
me14d...@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  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
> >
> > ah...@pitt.edu
> >
> > (724) 612-9940
> >
> >
> >
> > ----------
> > *From:* R-devel  on behalf of Durga
> Prasad
> > G me14d059 
> > *Sent:* Tuesday, August 1, 2023 4:18:20 AM
> > *To:* Martin Maechler ;
> r-devel@r-project.org
> > ; profjcn...@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 <
> > maech...@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
>

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

2023-08-16 Thread Lakshman, Aidan H
Hi Durga,

Just to add to previous comments--if you’re interested in the implementation of 
SVD in R, I’d recommend looking at the source code and/or sending general 
questions to the community slack, which may provide faster responses and help 
with double checking your work. Most of the perceived concerns you mention are 
related to the underlying LAPACK code--if you can verify that the code in 
LAPACK and/or R is incorrect, feel free to email r-devel with examples or 
discuss it on the slack.

The R implementation of SVD is available here, which essentially just wraps 
LAPACK code:
https://github.com/r-devel/r-svn/blob/b6394a83c0c12b12b6b1aceb05db0fd66227fd30/src/modules/lapack/Lapack.c#L111

The R method that invokes the C code is called here:
https://github.com/r-devel/r-svn/blob/b6394a83c0c12b12b6b1aceb05db0fd66227fd30/src/library/base/R/LAPACK.R#L19

The related LAPACK code is located here:
https://netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_gad8e0f1c83a78d3d4858eaaa88a1c5ab1.html

You can join the R-Contributors slack here:
https://contributor.r-project.org/slack

-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
ah...@pitt.edu
(724) 612-9940


From: Martin Maechler 
Date: Wednesday, August 16, 2023 at 04:40
To: Durga Prasad G me14d059 
Cc: Lakshman, Aidan H , Martin Maechler 
, r-devel@r-project.org , 
profjcn...@gmail.com 
Subject: Re: [Rd] Concerns with SVD -- and the Matrix Exponential
>>>>> Durga Prasad G me14d059
>>>>> on Wed, 16 Aug 2023 13:36:10 +0530 writes:

> Dear Martin, I am getting different responses from different officials of
> R-Software,

well, well, ..
Here on R-devel, we got two messages in addition to mine,
none by any "official" (even though John C Nash probably gets
the title of most senior professor still active on the R mailing list),
but both to the point,
notably Aidan Lakshman below showing you how you were unaware of
several things about the SVD, and confusing
positive semi definite matrices with arbitrary symmetric matrices.

> but those statements are contradicting with the statements
> discussed in your email.

I don't think so.  I think we all agreed (John Nash, me, Aidan Lakshman),
even though focussing on different aspects of your partly
incorrect claims.

Martin Maechler



> 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  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://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.ahl27.com%2F=05%7C01%7CAHL27%40pitt.edu%7C62ae1d193a224c3b1c3008db9e3469bb%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638277720201574436%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=gLTYfoHNU2p5JBF7ZHrUT%2BcFnZ%2BSJ6uqveFW3r6Sj84%3D=0<https://www.ahl27.com/>>
>>
>> Doctoral Fellow, Wright Lab 
<https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.wrightlabscience.com%2F=05%7C01%7CAHL27%40pitt.edu%7C62ae1d193a224c3b1c3008db9e3469bb%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638277720201574436%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=Wh7TaNdaRD%2B7euDlskDbzj7SBhPBt2y%2BFHcXkcYBM9M%3D=0<https://www.wrightlabscience.com/>>
>>
>> University of Pittsburgh School of Medicine
>>
>> Department of Biomedical Informatics
>>
>> ah...@pitt.edu
>>
>> (724) 612-9940
    >>
>>
>>
>> --
>> *From:* R-devel  on behalf of Durga Prasad
>> G me14d059 
>> *Sent:* Tuesday, August 1, 2023 4:18:20 AM

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

2023-08-16 Thread Durga Prasad G me14d059
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  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
>
> ah...@pitt.edu
>
> (724) 612-9940
>
>
>
> --
> *From:* R-devel  on behalf of Durga Prasad
> G me14d059 
> *Sent:* Tuesday, August 1, 2023 4:18:20 AM
> *To:* Martin Maechler ; r-devel@r-project.org
> ; profjcn...@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 <
> maech...@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=05%7C01%7Cahl27%40pitt.edu%7C8575b77db32345ca544b08db927ceae0%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638264837816871329%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=Y4mlFL%2FggLKd7FoIoY62esiFGUwukRG0YmELsJj7nd0%3D=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)
> > wr

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

2023-08-16 Thread Martin Maechler
>>>>> Durga Prasad G me14d059 
>>>>> on Wed, 16 Aug 2023 13:36:10 +0530 writes:

> Dear Martin, I am getting different responses from different officials of
> R-Software, 

well, well, ..
Here on R-devel, we got two messages in addition to mine,
none by any "official" (even though John C Nash probably gets
the title of most senior professor still active on the R mailing list),
but both to the point,
notably Aidan Lakshman below showing you how you were unaware of
several things about the SVD, and confusing
positive semi definite matrices with arbitrary symmetric matrices.

> but those statements are contradicting with the statements
> discussed in your email. 

I don't think so.  I think we all agreed (John Nash, me, Aidan Lakshman),
even though focussing on different aspects of your partly
incorrect claims.

Martin Maechler



> 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  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
>> 
>> ah...@pitt.edu
>> 
>> (724) 612-9940
>> 
>> 
    >> 
    >> ----------
>> *From:* R-devel  on behalf of Durga Prasad
>> G me14d059 
>> *Sent:* Tuesday, August 1, 2023 4:18:20 AM
>> *To:* Martin Maechler ; r-devel@r-project.org
>> ; profjcn...@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 <
>> maech...@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
>>

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

2023-08-01 Thread Lakshman, Aidan H
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

ah...@pitt.edu

(724) 612-9940




From: R-devel  on behalf of Durga Prasad G 
me14d059 
Sent: Tuesday, August 1, 2023 4:18:20 AM
To: Martin Maechler ; r-devel@r-project.org 
; profjcn...@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 
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=05%7C01%7Cahl27%40pitt.edu%7C8575b77db32345ca544b08db927ceae0%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638264837816871329%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=Y4mlFL%2FggLKd7FoIoY62esiFGUwukRG0YmELsJj7nd0%3D=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=05%7C01%7Cahl27%40pitt.edu%7C8575b77db32345ca544b08db927ceae0%7C9ef9f489e0a04eeb87cc3a526112fd0d%7C1%7C0%7C638264837817183809%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C=5%2FlssUGC6q7SUy0PY7gZCqWi0%2BXbNwZD0FaAgIcOWdY%3D=0<https://doi.org/10.1137/S00361445024180>
>
> It is BTW also cited on the Wikipedia page on the matrix
> exponential:
>
>
> ==> For this reas

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

2023-08-01 Thread Durga Prasad G me14d059
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 
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.
> >>
> >> 

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

2023-07-17 Thread 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.

__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel