R(matrix(1, 1, 1)) breaks - bug or desired?

[JanaJarecki]

Hi!

Is it desired that a one-dimensional matrix breaks the R() function? If so, maybe an error msg could help?

R(matrix(1, 1, 1))
> Fehler in U[, j] : Indizierung ausserhalb der Grenzen

Kind regards
Jana

[friendly]

Well, matrix(1,1,1) is a scalar, so R() should return 0 here.
matrix(1:3, 1,3) is a 1D matrix or vector, so R() should return 1.

Neither of these cases are presently trapped, so this is a reasonable TODO item

[john-d-fox]

Hi Michael, ----------------------------- John Fox, Professor Emeritus McMaster University Hamilton, Ontario, Canada Web: http::/socserv.mcmaster.ca/jfox

On May 3, 2020, at 8:22 PM, Michael Friendly ***@***.***> wrote: Well, matrix(1,1,1) is a scalar, so R() should return 0 here.

I don't agree: this is a 1 x 1 matrix with a nonzero element, not a scalar, so it should be of rank 1. It is an edge case, of course. Best, John

…

matrix(1:3, 1,3) is a 1D matrix or vector, so R() should return 1. Neither of these cases are presently trapped, so this is a reasonable TODO item — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[gmonette]

Hi, I think the problem is the classical dimension dropping problem in the QR function. Many of the matrix subsetting expressions need to have ‘,drop = FALSE’ added. I would offer but I’ve never touched the matlib code directly so I feel a bit shy! All the best, Georges PS. By the way qr(1)$rank returns 1. I.e. even if you feed it a scalar.

-- Georges Monette, PhD P.Stat. | Associate Professor, Faculty of Science, Department of Mathematics & Statistics | Office: Institute for Social Research: Dahdaleh Building 5055 | York University | 4700 Keele Street, Toronto, ON M3J 1P3 | E-Mail: georges@yorku.ca<mailto:georges@yorku.ca> Statistical Consulting Service Online Appointment Scheduler<http://www.appointmentquest.com/provider/2000199121> From: John Fox <notifications@github.com> Sent: May 3, 2020 9:01 PM To: friendly/matlib <matlib@noreply.github.com> Cc: Subscribed <subscribed@noreply.github.com> Subject: Re: [friendly/matlib] R(matrix(1, 1, 1)) breaks - bug or desired? (#34) Hi Michael, ----------------------------- John Fox, Professor Emeritus McMaster University Hamilton, Ontario, Canada Web: http::/socserv.mcmaster.ca/jfox

On May 3, 2020, at 8:22 PM, Michael Friendly ***@***.***> wrote: Well, matrix(1,1,1) is a scalar, so R() should return 0 here.

I don't agree: this is a 1 x 1 matrix with a nonzero element, not a scalar, so it should be of rank 1. It is an edge case, of course. Best, John

matrix(1:3, 1,3) is a 1D matrix or vector, so R() should return 1. Neither of these cases are presently trapped, so this is a reasonable TODO item — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#34 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAF76OSQQSOBQHNF34AHCMTRPYHU7ANCNFSM4MYKGRIA>.

[friendly]

> m <- matrix(1,1,1)
> qr(m)$rank
[1] 1
> qr(1:3)$rank
[1] 1

So, perhaps we should follow qr()

[john-d-fox]

Hi Michael, Since rank is the number of linearly independent rows or columns, this just seems correct -- that is not peculiar to qr(). Best, John

…

On May 3, 2020, at 11:01 PM, Michael Friendly ***@***.***> wrote: > m <- matrix(1,1,1) > qr(m)$rank [1] 1 > qr(1:3)$rank [1] 1 So, perhaps we should follow qr() — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[john-d-fox]

Hi Georges, Though you shouldn't hesitate to fix QR(), I can take a look at it sometime this week. There are other places where we had to set drop=FALSE to get correct behaviour for one-row or one-column cases. There are no true scalars in R, so the problem is technically a one-row or -column matrix (including, as here, a one-row and one-column matrix) becoming a vector. Best, John

…

On May 3, 2020, at 9:16 PM, Georges Monette ***@***.***> wrote: Hi, I think the problem is the classical dimension dropping problem in the QR function. Many of the matrix subsetting expressions need to have ‘,drop = FALSE’ added. I would offer but I’ve never touched the matlib code directly so I feel a bit shy! All the best, Georges PS. By the way qr(1)$rank returns 1. I.e. even if you feed it a scalar. -- Georges Monette, PhD P.Stat. | Associate Professor, Faculty of Science, Department of Mathematics & Statistics | Office: Institute for Social Research: Dahdaleh Building 5055 | York University | 4700 Keele Street, Toronto, ON M3J 1P3 | E-Mail: ***@***.******@***.***> Statistical Consulting Service Online Appointment Scheduler<http://www.appointmentquest.com/provider/2000199121> From: John Fox ***@***.***> Sent: May 3, 2020 9:01 PM To: friendly/matlib ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [friendly/matlib] R(matrix(1, 1, 1)) breaks - bug or desired? (#34) Hi Michael, ----------------------------- John Fox, Professor Emeritus McMaster University Hamilton, Ontario, Canada Web: http::/socserv.mcmaster.ca/jfox > On May 3, 2020, at 8:22 PM, Michael Friendly ***@***.***> wrote: > > > Well, matrix(1,1,1) is a scalar, so R() should return 0 here. I don't agree: this is a 1 x 1 matrix with a nonzero element, not a scalar, so it should be of rank 1. It is an edge case, of course. Best, John > matrix(1:3, 1,3) is a 1D matrix or vector, so R() should return 1. > > Neither of these cases are presently trapped, so this is a reasonable TODO item > > — > You are receiving this because you are subscribed to this thread. > Reply to this email directly, view it on GitHub, or unsubscribe. > — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#34 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAF76OSQQSOBQHNF34AHCMTRPYHU7ANCNFSM4MYKGRIA>. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[gmonette]

Hi John, Indeed, it never gets smaller than a vector in R. On a tangentially related issue, I wonder whether any of our packages will be affected by the new behaviour in version 4 of changing the default for ‘stringsAsFactors’ to FALSE. All the best, Georges From: John Fox <notifications@github.com> Sent: May 3, 2020 11:45 PM To: friendly/matlib <matlib@noreply.github.com> Cc: Georges A Monette <georges@yorku.ca>; Comment <comment@noreply.github.com> Subject: Re: [friendly/matlib] R(matrix(1, 1, 1)) breaks - bug or desired? (#34) Hi Georges, Though you shouldn't hesitate to fix QR(), I can take a look at it sometime this week. There are other places where we had to set drop=FALSE to get correct behaviour for one-row or one-column cases. There are no true scalars in R, so the problem is technically a one-row or -column matrix (including, as here, a one-row and one-column matrix) becoming a vector. Best, John

On May 3, 2020, at 9:16 PM, Georges Monette ***@***.***> wrote: Hi, I think the problem is the classical dimension dropping problem in the QR function. Many of the matrix subsetting expressions need to have ‘,drop = FALSE’ added. I would offer but I’ve never touched the matlib code directly so I feel a bit shy! All the best, Georges PS. By the way qr(1)$rank returns 1. I.e. even if you feed it a scalar. -- Georges Monette, PhD P.Stat. | Associate Professor, Faculty of Science, Department of Mathematics & Statistics | Office: Institute for Social Research: Dahdaleh Building 5055 | York University | 4700 Keele Street, Toronto, ON M3J 1P3 | E-Mail: ***@***.******@***.***> Statistical Consulting Service Online Appointment Scheduler<http://www.appointmentquest.com/provider/2000199121> From: John Fox ***@***.***> Sent: May 3, 2020 9:01 PM To: friendly/matlib ***@***.***> Cc: Subscribed ***@***.***> Subject: Re: [friendly/matlib] R(matrix(1, 1, 1)) breaks - bug or desired? (#34) Hi Michael, ----------------------------- John Fox, Professor Emeritus McMaster University Hamilton, Ontario, Canada Web: http::/socserv.mcmaster.ca/jfox > On May 3, 2020, at 8:22 PM, Michael Friendly ***@***.***> wrote: > > > Well, matrix(1,1,1) is a scalar, so R() should return 0 here. I don't agree: this is a 1 x 1 matrix with a nonzero element, not a scalar, so it should be of rank 1. It is an edge case, of course. Best, John > matrix(1:3, 1,3) is a 1D matrix or vector, so R() should return 1. > > Neither of these cases are presently trapped, so this is a reasonable TODO item > > — > You are receiving this because you are subscribed to this thread. > Reply to this email directly, view it on GitHub, or unsubscribe. > — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#34 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAF76OSQQSOBQHNF34AHCMTRPYHU7ANCNFSM4MYKGRIA>. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub<#34 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAF76OT3QSVBH63LTYZO2ZTRPY223ANCNFSM4MYKGRIA>.

[john-d-fox]

Hi Georges, It's slightly off topic, but if by "our" you mean packages maintained not necessarily by all of us, then the answer is yes. There are small implications for both the car and Rcmdr packages. It's possible that others are affected as well, but I'm unaware of them. It's certainly worth checking. Best, John

…

On May 4, 2020, at 9:48 AM, Georges Monette ***@***.***> wrote: Hi John, Indeed, it never gets smaller than a vector in R. On a tangentially related issue, I wonder whether any of our packages will be affected by the new behaviour in version 4 of changing the default for ‘stringsAsFactors’ to FALSE. All the best, Georges From: John Fox ***@***.***> Sent: May 3, 2020 11:45 PM To: friendly/matlib ***@***.***> Cc: Georges A Monette ***@***.***>; Comment ***@***.***> Subject: Re: [friendly/matlib] R(matrix(1, 1, 1)) breaks - bug or desired? (#34) Hi Georges, Though you shouldn't hesitate to fix QR(), I can take a look at it sometime this week. There are other places where we had to set drop=FALSE to get correct behaviour for one-row or one-column cases. There are no true scalars in R, so the problem is technically a one-row or -column matrix (including, as here, a one-row and one-column matrix) becoming a vector. Best, John > On May 3, 2020, at 9:16 PM, Georges Monette ***@***.***> wrote: > > > Hi, > I think the problem is the classical dimension dropping problem in the QR function. Many of the matrix subsetting expressions need to have ‘,drop = FALSE’ added. I would offer but I’ve never touched the matlib code directly so I feel a bit shy! > All the best, > Georges > > PS. By the way qr(1)$rank returns 1. I.e. even if you feed it a scalar. > > -- > Georges Monette, PhD P.Stat. | Associate Professor, Faculty of Science, Department of Mathematics & Statistics | Office: Institute for Social Research: Dahdaleh Building 5055 | York University | 4700 Keele Street, Toronto, ON M3J 1P3 | E-Mail: ***@***.******@***.***> > > Statistical Consulting Service Online Appointment Scheduler<http://www.appointmentquest.com/provider/2000199121> > > From: John Fox ***@***.***> > Sent: May 3, 2020 9:01 PM > To: friendly/matlib ***@***.***> > Cc: Subscribed ***@***.***> > Subject: Re: [friendly/matlib] R(matrix(1, 1, 1)) breaks - bug or desired? (#34) > > Hi Michael, > ----------------------------- > John Fox, Professor Emeritus > McMaster University > Hamilton, Ontario, Canada > Web: http::/socserv.mcmaster.ca/jfox > > > On May 3, 2020, at 8:22 PM, Michael Friendly ***@***.***> wrote: > > > > > > Well, matrix(1,1,1) is a scalar, so R() should return 0 here. > > I don't agree: this is a 1 x 1 matrix with a nonzero element, not a scalar, so it should be of rank 1. It is an edge case, of course. > > Best, > John > > > matrix(1:3, 1,3) is a 1D matrix or vector, so R() should return 1. > > > > Neither of these cases are presently trapped, so this is a reasonable TODO item > > > > — > > You are receiving this because you are subscribed to this thread. > > Reply to this email directly, view it on GitHub, or unsubscribe. > > > > > > — > You are receiving this because you are subscribed to this thread. > Reply to this email directly, view it on GitHub<#34 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAF76OSQQSOBQHNF34AHCMTRPYHU7ANCNFSM4MYKGRIA>. > — > You are receiving this because you commented. > Reply to this email directly, view it on GitHub, or unsubscribe. > — You are receiving this because you commented. Reply to this email directly, view it on GitHub<#34 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAF76OT3QSVBH63LTYZO2ZTRPY223ANCNFSM4MYKGRIA>. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[friendly]

R(X) is just short for QR(X)$rank. And recall that QR() was designed to make the computation visible to inspection, both at the code and computation levels.

A simple fix to QR() avoids this problem of R() generating an error when ncol{X} < 2.

  if (ncol(U)>1) {                 # trap potential error
    for (j in 2:ncol(U)){
      for (k in 1:(j - 1)){
        U[, j] <- U[, j] - as.vector(X[, j] %*% E[, k]) * E[, k]
      }
      len.U.j <- length(U[, j])
      if (len.U.j > tol) E[, j] <- U[, j]/len.U.j
    }
  }

But I don't think this is a general solution, at least if we want QR to match qr in such edge cases.

> QR(matrix(1:3,1,3))
$Q
     [,1] [,2] [,3]
[1,]   -1    0    0

$R
     [,1] [,2] [,3]
[1,]   -1   -2   -3
[2,]    0    0    0
[3,]    0    0    0

$rank
[1] 1

> qr(matrix(1:3,1,3))
$qr
     [,1] [,2] [,3]
[1,]    1    2    3

$rank
[1] 1

[john-d-fox]

Hi Michael, If I remember right, I wrote the initial version of QR() and Georges improved it. Unless Georges wants to do it, I can take a look at this, probably sometime this week. Best, John

…

On May 5, 2020, at 10:27 AM, Michael Friendly ***@***.***> wrote: R(X) is just short for QR(X)$rank. And recall that QR() was designed to make the computation visible to inspection, both at the code and computation levels. A simple fix to QR() avoids this problem of R() generating an error when ncol{X} < 2. if (ncol(U)>1) { # trap potential error for (j in 2:ncol(U)){ for (k in 1:(j - 1)){ U[, j] <- U[, j] - as.vector(X[, j] %*% E[, k]) * E[, k] } len.U.j <- length(U[, j]) if (len.U.j > tol) E[, j] <- U[, j]/len.U.j } } But I don't think this is a general solution, at least if we want QR to match qr in such edge cases. > QR(matrix(1:3,1,3)) $Q [,1] [,2] [,3] [1,] -1 0 0 $R [,1] [,2] [,3] [1,] -1 -2 -3 [2,] 0 0 0 [3,] 0 0 0 $rank [1] 1 > qr(matrix(1:3,1,3)) $qr [,1] [,2] [,3] [1,] 1 2 3 $rank [1] 1 — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[gmonette]

Hi, I can give it a go this week since I’m working the problem of correctly identifying the rank of a spline for gspline, there will be synergies! All the best, Georges

-- Georges Monette, PhD P.Stat. | Associate Professor, Faculty of Science, Department of Mathematics & Statistics | Office: Institute for Social Research: Dahdaleh Building 5055 | York University | 4700 Keele Street, Toronto, ON M3J 1P3 | E-Mail: georges@yorku.ca<mailto:georges@yorku.ca> Statistical Consulting Service Online Appointment Scheduler<http://www.appointmentquest.com/provider/2000199121> From: John Fox <notifications@github.com> Sent: May 5, 2020 11:17 AM To: friendly/matlib <matlib@noreply.github.com> Cc: Georges A Monette <georges@yorku.ca>; Comment <comment@noreply.github.com> Subject: Re: [friendly/matlib] R(matrix(1, 1, 1)) breaks - bug or desired? (#34) Hi Michael, If I remember right, I wrote the initial version of QR() and Georges improved it. Unless Georges wants to do it, I can take a look at this, probably sometime this week. Best, John

On May 5, 2020, at 10:27 AM, Michael Friendly ***@***.***> wrote: R(X) is just short for QR(X)$rank. And recall that QR() was designed to make the computation visible to inspection, both at the code and computation levels. A simple fix to QR() avoids this problem of R() generating an error when ncol{X} < 2. if (ncol(U)>1) { # trap potential error for (j in 2:ncol(U)){ for (k in 1:(j - 1)){ U[, j] <- U[, j] - as.vector(X[, j] %*% E[, k]) * E[, k] } len.U.j <- length(U[, j]) if (len.U.j > tol) E[, j] <- U[, j]/len.U.j } } But I don't think this is a general solution, at least if we want QR to match qr in such edge cases. > QR(matrix(1:3,1,3)) $Q [,1] [,2] [,3] [1,] -1 0 0 $R [,1] [,2] [,3] [1,] -1 -2 -3 [2,] 0 0 0 [3,] 0 0 0 $rank [1] 1 > qr(matrix(1:3,1,3)) $qr [,1] [,2] [,3] [1,] 1 2 3 $rank [1] 1 — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub<#34 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAF76OVO6R4QJRYVWGFOLFLRQAUYTANCNFSM4MYKGRIA>.

[john-d-fox]

Hi, I took a look at this before seeing Georges's message, and I don't see what the problem is with Michael's fix. That is, the differences with qr() are inessential (arbitrary compensating sign changes and inclusion/suppression of 0 rows). Here are some examples, including three edge cases (case B is Michael's):

A <- matrix(1:9, 3, 3) B <- matrix(1:3, 1) C <- matrix(1:3, ncol=1) D <- matrix(10, 1, 1) A

[,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9

(QRA <- QR(A))

$Q [,1] [,2] [,3] [1,] -0.2672612 -0.8728716 0 [2,] -0.5345225 -0.2182179 0 [3,] -0.8017837 0.4364358 0 $R [,1] [,2] [,3] [1,] -3.741657 -8.552360 -13.363062 [2,] 0.000000 -1.963961 -3.927922 [3,] 0.000000 0.000000 0.000000 $rank [1] 2

QRA$Q %*% QRA$R

[,1] [,2] [,3] [1,] 1 4 7 [2,] 2 5 8 [3,] 3 6 9

qrA <- qr(A) qr.Q(qrA)

[,1] [,2] [,3] [1,] -0.2672612 0.8728716 0.4082483 [2,] -0.5345225 0.2182179 -0.8164966 [3,] -0.8017837 -0.4364358 0.4082483

qr.R(qrA)

[,1] [,2] [,3] [1,] -3.741657 -8.552360 -1.336306e+01 [2,] 0.000000 1.963961 3.927922e+00 [3,] 0.000000 0.000000 1.776357e-15

B

[,1] [,2] [,3] [1,] 1 2 3

(QRB <- QR(B))

$Q [,1] [,2] [,3] [1,] -1 0 0 $R [,1] [,2] [,3] [1,] -1 -2 -3 [2,] 0 0 0 [3,] 0 0 0 $rank [1] 1

QRB$Q %*% QRB$R

[,1] [,2] [,3] [1,] 1 2 3

qrB <- qr(B) qr.Q(qrB)

[,1] [1,] 1

qr.R(qrB)

[,1] [,2] [,3] [1,] 1 2 3

C

[,1] [1,] 1 [2,] 2 [3,] 3

(QRC <- QR(C))

$Q [,1] [1,] -0.2672612 [2,] -0.5345225 [3,] -0.8017837 $R [,1] [1,] -3.741657 $rank [1] 1

QRC$Q %*% QRC$R

[,1] [1,] 1 [2,] 2 [3,] 3

qrC <- qr(C) qr.Q(qrC)

[,1] [1,] -0.2672612 [2,] -0.5345225 [3,] -0.8017837

qr.R(qrC)

[,1] [1,] -3.741657

D

[,1] [1,] 10

(QRD <- QR(D))

$Q [,1] [1,] -1 $R [,1] [1,] -10 $rank [1] 1

QRD$Q %*% QRD$R

[,1] [1,] 10

qrD <- qr(D) qr.Q(qrD)

[,1] [1,] 1

qr.R(qrD)

[,1] [1,] 10 All this looks right to me. Best, John

…

On May 5, 2020, at 3:58 PM, Georges Monette ***@***.***> wrote: Hi, I can give it a go this week since I’m working the problem of correctly identifying the rank of a spline for gspline, there will be synergies! All the best, Georges -- Georges Monette, PhD P.Stat. | Associate Professor, Faculty of Science, Department of Mathematics & Statistics | Office: Institute for Social Research: Dahdaleh Building 5055 | York University | 4700 Keele Street, Toronto, ON M3J 1P3 | E-Mail: ***@***.******@***.***> Statistical Consulting Service Online Appointment Scheduler<http://www.appointmentquest.com/provider/2000199121> From: John Fox ***@***.***> Sent: May 5, 2020 11:17 AM To: friendly/matlib ***@***.***> Cc: Georges A Monette ***@***.***>; Comment ***@***.***> Subject: Re: [friendly/matlib] R(matrix(1, 1, 1)) breaks - bug or desired? (#34) Hi Michael, If I remember right, I wrote the initial version of QR() and Georges improved it. Unless Georges wants to do it, I can take a look at this, probably sometime this week. Best, John > On May 5, 2020, at 10:27 AM, Michael Friendly ***@***.***> wrote: > > > R(X) is just short for QR(X)$rank. And recall that QR() was designed to make the computation visible to inspection, both at the code and computation levels. > > A simple fix to QR() avoids this problem of R() generating an error when ncol{X} < 2. > > if (ncol(U)>1) { # trap potential error > for (j in 2:ncol(U)){ > for (k in 1:(j - 1)){ > U[, j] <- U[, j] - as.vector(X[, j] %*% E[, k]) * E[, k] > } > len.U.j <- length(U[, j]) > if (len.U.j > tol) E[, j] <- U[, j]/len.U.j > } > } > > But I don't think this is a general solution, at least if we want QR to match qr in such edge cases. > > > QR(matrix(1:3,1,3)) > $Q > [,1] [,2] [,3] > [1,] -1 0 0 > > $R > [,1] [,2] [,3] > [1,] -1 -2 -3 > [2,] 0 0 0 > [3,] 0 0 0 > > $rank > [1] 1 > > > qr(matrix(1:3,1,3)) > $qr > [,1] [,2] [,3] > [1,] 1 2 3 > > $rank > [1] 1 > > — > You are receiving this because you commented. > Reply to this email directly, view it on GitHub, or unsubscribe. > — You are receiving this because you commented. Reply to this email directly, view it on GitHub<#34 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AAF76OVO6R4QJRYVWGFOLFLRQAUYTANCNFSM4MYKGRIA>. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.