search

friendly / matlib

Matrix Functions for Teaching and Learning Linear Algebra and Multivariate Statistics
http://friendly.github.io/matlib/
65 stars 16 forks source link

Consider adding arguments cycles and max.denominator to formatNumbers() #42

Closed TengMCing closed 1 year ago

TengMCing commented 3 years ago

Hi,

when I am using 

A <- matrix(c(1/5*12^4, 0, 48, 0, 1,
              864, 108, 6, 1, 0,
              -864, 108, -6, 1, 0,
              0, 144, 0, 1, 0),
            byrow = TRUE,
            nrow = 4)

b <- matrix(c(1,0,0,0), nrow = 4)

matlib::gaussianElimination(A, b, verbose = TRUE, fractions = TRUE)

the output is something like

Initial matrix:
     [,1]    [,2]    [,3]    [,4]    [,5]    [,6]   
[1,] 20736/5       0      48       0       1       1
[2,]     864     108       6       1       0       0
[3,]    -864     108      -6       1       0       0
[4,]       0     144       0       1       0       0

row: 1 

 multiply row 1 by 0 
     [,1]  [,2]  [,3]  [,4]  [,5]  [,6] 
[1,]     1     0 5/432     0     0     0
[2,]   864   108     6     1     0     0
[3,]  -864   108    -6     1     0     0
[4,]     0   144     0     1     0     0

.
.
.
 multiply row 4 by 1/16 and subtract from row 3 
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    0    0    0    0
[2,]    0    1    0    0    0    0
[3,]    0    0    1    0 5/96 5/96
[4,]    0    0    0    1    0    0

which is incorrect. The first step round the value in row 1 col 5 and 6 to zero. However, if I set fractions=FALSE, it gives the correct answer. So it has to be something wrong with the fractions.

After inspecting the source code, I find there is a potential improvement that can be made in util.R:

formatNumbers <- function(x, fractions=FALSE, tol=sqrt(.Machine$double.eps)){
    ret <- if (fractions) MASS::fractions(x) else round(x, round(abs(log(tol, 10))))
    ret
}

Using the default values of MASS::fractions sometimes could be insufficient. I think you could allow users to pass arguments cycles and max.determinator to MASS::fractions, or set the default values a bit larger to overcome the issue.

john-d-fox commented 3 years ago

Hi Patrick,

Thanks for pointing out this problem. I don't see an easy fix. MASS::fractions() is called at a number of different points in the package code.

We could pass through a ... argument, since fractions() has one and thus irrelevant specific arguments would then be (properly) ignored, but typically fractions() isn't called directly from a top-level function in matlib.

Another solution would be to introduce something like a fractions.precision option(), but I'd rather stick with functional programming.

Finally, we could introduce out own Fractions() function, which calls MASS::fractions with different defaults for cycles and max.denominator, but the values we pick would still be hard-coded and might also fail.

Of these possibilities (can you think of others?), maybe the last is best.

Thanks again, John

John Fox, Professor Emeritus McMaster University Hamilton, Ontario, Canada web: https://socialsciences.mcmaster.ca/jfox/

On 2021-08-29 3:02 a.m., Weihao (Patrick) Li wrote:

Hi,

when I am using

A <- matrix(c(1/5*12^4,0,48,0,1, 864,108,6,1,0, -864,108,-6,1,0, 0,144,0,1,0), byrow = TRUE, nrow = 4)

b <- matrix(c(1,0,0,0),nrow = 4)

matlib::gaussianElimination(A,b,verbose = TRUE,fractions = TRUE)

the output is something like

|Initial matrix: [,1] [,2] [,3] [,4] [,5] [,6] [1,] 20736/5 0 48 0 1 1 [2,] 864 108 6 1 0 0 [3,] -864 108 -6 1 0 0 [4,] 0 144 0 1 0 0 row: 1 multiply row 1 by 0 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 0 0 [2,] 864 108 6 1 0 0 [3,] -864 108 -6 1 0 0 [4,] 0 144 0 1 0 0 .. . .. multiply row 4 by 1/16 and subtract from row 3 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 0 0 0 0 [2,] 0 1 0 0 0 0 [3,] 0 0 1 0 5/96 5/96 [4,] 0 0 0 1 0 0 |

which is incorrect. The first step round the value in row 1 col 5 and 6 to zero. However, if I set |fractions=FALSE|, it gives the correct answer. So it has to be something wrong with the fractions.

After inspecting the source code, I find there is a potential improvement that can be made in util.R:

|formatNumbers <- function(x, fractions=FALSE, tol=sqrt(.Machine$double.eps)){ ret <- if (fractions) MASS::fractions(x) else round(x, round(abs(log(tol, 10)))) ret } |

Using the default values of |MASS::fractions| sometimes could be insufficient. I think you could allow users to pass arguments |cycles| and |max.determinator| to |MASS::fractions|, or set the default values a bit larger to overcome the issue.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/friendly/matlib/issues/42, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLSAQR7UOP6CPIQFXJVA63T7HLRHANCNFSM5C75XVLA. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

TengMCing commented 3 years ago

I guess a quick fix will be to introduce a function factory say 

Fractions <- function(cycles = 10, max.denominator = 2000) {
  customFractions <- TRUE
  attr(customFractions, "formatFunc") <- function(x) {
    MASS::fractions(x, cycles = cycles, max.denominator)
  }
  customFractions
}

and the formatNumbers() could be modified to

formatNumbers <- function(x, fractions=FALSE, tol=sqrt(.Machine$double.eps)){
  ret <- if (fractions) {
    if ("formatFunc" %in% names(attributes(fractions))) {
      attr(fractions, "formatFunc")(x)
    } else {
      MASS::fractions(x)
    }
  } else round(x, round(abs(log(tol, 10))))
  ret
}

then allow the user to set fractions =TRUE, fractions = FALSE, or define its own fraction function via the function factory fractions = Fractions(cycles = mycycles, max.denominator = mymax.denominator).

This will not affect the main functionality of the package since the returned object of Fractions() is still a logical value, but it allows the user to control the fraction formatting.

john-d-fox commented 3 years ago

Hi,

That's clever and should work (except that MASS::fractions() isn't always called via formatNumbers()). Also, documenting the solution in an understandable manner for users would probably be challenging, and perhaps is too much for what's likely an unusual problem.

I'm inclined simply to replace MASS::fractions with Fractions(), and despite aversion to a non-functional solution, to allow a fractions option() that can set cycles and max.denominator.

Best, John On 2021-08-29 9:22 a.m., Weihao (Patrick) Li wrote:

I guess a quick fix will be to introduce a function factory say

|Fractions <- function(cycles = 10, max.denominator = 2000) { customFractions <- TRUE attr(customFractions, "formatFunc") <- function(x = NULL) { MASS::fractions(x, cycles = cycles, max.denominator) } customFractions } |

and the |formatNumbers()| could be modified to

|formatNumbers <- function(x, fractions=FALSE, tol=sqrt(.Machine$double.eps)){ ret <- if (fractions) { if ("formatFunc" %in% names(attributes(fractions))) { attr(fractions, "formatFunc")(x) } else { MASS::fractions(x) } } else round(x, round(abs(log(tol, 10)))) ret } |

then allow the user to set |fractions =TRUE|, |fractions = FALSE|, or define its own fraction function via the function factory |fractions = Fractions(cycles = mycycles, max.denominator = mymax.denominator)|.

This will not affect the main functionality of the package since the returned object of |Fractions()| is still a logical value, but it allows the user to control the fraction formatting.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/friendly/matlib/issues/42#issuecomment-907790741, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLSAQROLSPZ55D7BDCJ4W3T7IX73ANCNFSM5C75XVLA. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

john-d-fox commented 3 years ago

I've now implemented but not committed code that does the following:

------ snip --------

options(fractions=list(cycles=1000, max.denominator=10^6)) # can specify both or neither gaussianElimination(A, b, verbose = TRUE, fractions = TRUE)

Initial matrix: [,1] [,2] [,3] [,4] [,5] [,6] [1,] 20736/5 0 48 0 1 1 [2,] 864 108 6 1 0 0 [3,] -864 108 -6 1 0 0 [4,] 0 144 0 1 0 0

row: 1

multiply row 1 by 5/20736 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 5/20736 5/20736 [2,] 864 108 6 1 0 0 [3,] -864 108 -6 1 0 0 [4,] 0 144 0 1 0 0

multiply row 1 by 864 and subtract from row 2 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 5/20736 5/20736 [2,] 0 108 -4 1 -5/24 -5/24 [3,] -864 108 -6 1 0 0 [4,] 0 144 0 1 0 0

multiply row 1 by 864 and add to row 3 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 5/20736 5/20736 [2,] 0 108 -4 1 -5/24 -5/24 [3,] 0 108 4 1 5/24 5/24 [4,] 0 144 0 1 0 0

row: 2

exchange rows 2 and 4 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 5/20736 5/20736 [2,] 0 144 0 1 0 0 [3,] 0 108 4 1 5/24 5/24 [4,] 0 108 -4 1 -5/24 -5/24

multiply row 2 by 1/144 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 5/20736 5/20736 [2,] 0 1 0 1/144 0 0 [3,] 0 108 4 1 5/24 5/24 [4,] 0 108 -4 1 -5/24 -5/24

multiply row 2 by 108 and subtract from row 3 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 5/20736 5/20736 [2,] 0 1 0 1/144 0 0 [3,] 0 0 4 1/4 5/24 5/24 [4,] 0 108 -4 1 -5/24 -5/24

multiply row 2 by 108 and subtract from row 4 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 5/20736 5/20736 [2,] 0 1 0 1/144 0 0 [3,] 0 0 4 1/4 5/24 5/24 [4,] 0 0 -4 1/4 -5/24 -5/24

row: 3

multiply row 3 by 1/4 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 5/432 0 5/20736 5/20736 [2,] 0 1 0 1/144 0 0 [3,] 0 0 1 1/16 5/96 5/96 [4,] 0 0 -4 1/4 -5/24 -5/24

multiply row 3 by 5/432 and subtract from row 1 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 0 -5/6912 -5/13824 -5/13824 [2,] 0 1 0 1/144 0 0 [3,] 0 0 1 1/16 5/96 5/96 [4,] 0 0 -4 1/4 -5/24 -5/24

multiply row 3 by 4 and add to row 4 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 0 -5/6912 -5/13824 -5/13824 [2,] 0 1 0 1/144 0 0 [3,] 0 0 1 1/16 5/96 5/96 [4,] 0 0 0 1/2 0 0

row: 4

multiply row 4 by 2 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 0 -5/6912 -5/13824 -5/13824 [2,] 0 1 0 1/144 0 0 [3,] 0 0 1 1/16 5/96 5/96 [4,] 0 0 0 1 0 0

multiply row 4 by 5/6912 and add to row 1 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 0 0 -5/13824 -5/13824 [2,] 0 1 0 1/144 0 0 [3,] 0 0 1 1/16 5/96 5/96 [4,] 0 0 0 1 0 0

multiply row 4 by 1/144 and subtract from row 2 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 0 0 -5/13824 -5/13824 [2,] 0 1 0 0 0 0 [3,] 0 0 1 1/16 5/96 5/96 [4,] 0 0 0 1 0 0

multiply row 4 by 1/16 and subtract from row 3 [,1] [,2] [,3] [,4] [,5] [,6] [1,] 1 0 0 0 -5/13824 -5/13824 [2,] 0 1 0 0 0 0 [3,] 0 0 1 0 5/96 5/96 [4,] 0 0 0 1 0 0

------ snip --------

Before committing and documenting this solution, thoughts from anyone?

John

On 2021-08-29 7:18 p.m., John Fox wrote:

Hi,

That's clever and should work (except that MASS::fractions() isn't always called via formatNumbers()). Also, documenting the solution in an understandable manner for users would probably be challenging, and perhaps is too much for what's likely an unusual problem.

I'm inclined simply to replace MASS::fractions with Fractions(), and despite aversion to a non-functional solution, to allow a fractions option() that can set cycles and max.denominator.

Best,  John On 2021-08-29 9:22 a.m., Weihao (Patrick) Li wrote:

I guess a quick fix will be to introduce a function factory say

|Fractions <- function(cycles = 10, max.denominator = 2000) { customFractions <- TRUE attr(customFractions, "formatFunc") <- function(x = NULL) { MASS::fractions(x, cycles = cycles, max.denominator) } customFractions } |

and the |formatNumbers()| could be modified to

|formatNumbers <- function(x, fractions=FALSE, tol=sqrt(.Machine$double.eps)){ ret <- if (fractions) { if ("formatFunc" %in% names(attributes(fractions))) { attr(fractions, "formatFunc")(x) } else { MASS::fractions(x) } } else round(x, round(abs(log(tol, 10)))) ret } |

then allow the user to set |fractions =TRUE|, |fractions = FALSE|, or define its own fraction function via the function factory |fractions = Fractions(cycles = mycycles, max.denominator = mymax.denominator)|.

This will not affect the main functionality of the package since the returned object of |Fractions()| is still a logical value, but it allows the user to control the fraction formatting.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/friendly/matlib/issues/42#issuecomment-907790741, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLSAQROLSPZ55D7BDCJ4W3T7IX73ANCNFSM5C75XVLA.

Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

TengMCing commented 3 years ago

Looks really good! I will close this issue as it has been resolved.

john-d-fox commented 3 years ago

On Aug 30, 2021, at 11:44 PM, Weihao (Patrick) Li @.***> wrote:



Looks really good! I will close this issue as it has been resolved.

Not quite yet please. I’d like to hear what Michael and Phil have to say. 

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android.

friendly commented 3 years ago

@john-d-fox Given that this problem seems to be relatively rare, your solution seems more than adequate. I recommend going ahead with this commit and adding some documentation on the options(fractions=).

Perhaps @philchalmers could also look this over

philchalmers commented 3 years ago

The proposed solution seems reasonable to me, though the use of options() in packages always feels a little awkward.

As a follow up, is there any reason not to increase the values of cycles and max.denomom beyond the MASS::fractions() values by default? The fraction calls are only used in verbose mode in the console, in which case users will be less interested in computational performance anyway.

Should we create a new temporary branch to test John's code prior to merging it into the master branch?

john-d-fox commented 3 years ago

Hi, On 2021-08-31 10:03 a.m., Phil Chalmers wrote:

The proposed solution seems reasonable to me, though the use of |options()| in packages always feels a little awkward.

Agreed!

As a follow up, is there any reason not to increase the values of |cycles| and |max.denomom| beyond the |MASS::fractions()| values by default?

Venables and Ripley are pretty smart guys and I bet they put some thought into the defaults, but the way I wrote Fractions() makes it easy to change the defaults.

The fraction calls are only used in verbose mode in the console, in which case users will be less interested in computational performance anyway.

Should we create a new temporary branch to test John's code prior to merging it into the |master| branch?

I'm having trouble pushing my changes. Apparently, I used a password with the matlib github repository, even though I have an RSA key set up. and passwords are now defunct.

Anyway, the changes seem to work fine in my local copy and are now also in the roxygen documentation.

John

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/friendly/matlib/issues/42#issuecomment-909268583, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLSAQUZSH2RUW2KGAVZUGLT7TOMHANCNFSM5C75XVLA. Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.

friendly commented 1 year ago

Looking over this old thread, I'm inclined to merge this into master and close this issue...

It passes R CMD check, but I wish we had some unit testing here...