friendly / matlib

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

latexMatrix arithmetic: trivial simplifications? #58

Closed friendly closed 2 months ago

friendly commented 2 months ago

In my vignette example of linear hypotheses, the following generates the C & B matrices:

C <- latexMatrix(matrix(c(0, 1, 0, 0,
                          0, 0, 1, 0), nrow=2, byrow=TRUE), 
                 matrix = "bmatrix")
B <- latexMatrix('\\beta', ncol = 3, nrow=4, 
                 comma=TRUE, prefix.col = 'y_',
                 zero.based=c(TRUE, FALSE))
C %*%B

which generates a correct result, but one that is hard to look at:

\begin{bmatrix}  
0 \cdot \beta_{0,y_{1}} + 1 \cdot \beta_{1,y_{1}} + 0 \cdot \beta_{2,y_{1}} + 0 \cdot \beta_{3,y_{1}} & 0 \cdot \beta_{0,y_{2}} + 1 \cdot \beta_{1,y_{2}} + 0 \cdot \beta_{2,y_{2}} + 0 \cdot \beta_{3,y_{2}} & 0 \cdot \beta_{0,y_{3}} + 1 \cdot \beta_{1,y_{3}} + 0 \cdot \beta_{2,y_{3}} + 0 \cdot \beta_{3,y_{3}} \\ 
0 \cdot \beta_{0,y_{1}} + 0 \cdot \beta_{1,y_{1}} + 1 \cdot \beta_{2,y_{1}} + 0 \cdot \beta_{3,y_{1}} & 0 \cdot \beta_{0,y_{2}} + 0 \cdot \beta_{1,y_{2}} + 1 \cdot \beta_{2,y_{2}} + 0 \cdot \beta_{3,y_{2}} & 0 \cdot \beta_{0,y_{3}} + 0 \cdot \beta_{1,y_{3}} + 1 \cdot \beta_{2,y_{3}} + 0 \cdot \beta_{3,y_{3}} \\ 
\end{bmatrix}

and looks like this:

image

Trivial simplification of this could:

Just removing the 0 * terms:

\begin{bmatrix}  
  1 \cdot \beta_{1,y_{1}} &  1 \cdot \beta_{1,y_{2}}  &  1 \cdot \beta_{1,y_{3}} \\ 
  1 \cdot \beta_{2,y_{1}} &  1 \cdot \beta_{2,y_{2}}  &  1 \cdot \beta_{2,y_{3}} \\ 
\end{bmatrix}

Is this possible? An operator function like %*% can't take other arguments, but could it handle an option in the environment?

philchalmers commented 2 months ago

It's certainly possible to do this. I did something similar in dev/Arith2.R for the kronecker product, and I was glad to see this idea made into John's implementation as well. Perhaps the best way to do this is to add a print(..., sparse = TRUE) to detect any zero multiplications in defined matrices and have them drop out?

john-d-fox commented 2 months ago

I think that sparse = TRUE is different: replacing 0s with blanks. That's a worthwhile option, but it's different from, and simpler than, replacing expressions that evaluate to 0 with 0 (which the kronecker() method also does) and removing multiplication by 1.

I made good progress on a simplify() function yesterday evening. I have appointments today, so probably won't get much done, but should have a preliminary implementation soon.

simplify() could optionally be called by print.latexMatrix() or used directly.

John

On 2024-08-21 11:07 p.m., Phil Chalmers wrote:

Caution: External email.

It's certainly possible to do this. I did something similar in |dev/ Arith2.R| for the kronecker product, and I was glad to see this idea made into John's implementation as well. Perhaps the best way to do this is to add a |print(..., sparse = TRUE)| to detect any zero multiplications in defined matrices and have them drop out?

— Reply to this email directly, view it on GitHub <https://github.com/ friendly/matlib/issues/58#issuecomment-2303586612>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ ADLSAQU67XI3LYXXCNX4KBTZSVIWFAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBTGU4DMNRRGI>. You are receiving this because you are subscribed to this thread.Message ID: @.***>

friendly commented 2 months ago

Another option is to have a parallel function, matmult(A, B), that takes a simplify argument More on this later

Get Outlook for iOShttps://aka.ms/o0ukef


From: John Fox @.> Sent: Thursday, August 22, 2024 1:28:29 PM To: friendly/matlib @.> Cc: Michael L Friendly @.>; Author @.> Subject: Re: [friendly/matlib] latexMatrix arithmetic: trivial simplifications? (Issue #58)

I think that sparse = TRUE is different: replacing 0s with blanks. That's a worthwhile option, but it's different from, and simpler than, replacing expressions that evaluate to 0 with 0 (which the kronecker() method also does) and removing multiplication by 1.

I made good progress on a simplify() function yesterday evening. I have appointments today, so probably won't get much done, but should have a preliminary implementation soon.

simplify() could optionally be called by print.latexMatrix() or used directly.

John

On 2024-08-21 11:07 p.m., Phil Chalmers wrote:

Caution: External email.

It's certainly possible to do this. I did something similar in |dev/ Arith2.R| for the kronecker product, and I was glad to see this idea made into John's implementation as well. Perhaps the best way to do this is to add a |print(..., sparse = TRUE)| to detect any zero multiplications in defined matrices and have them drop out?

— Reply to this email directly, view it on GitHub https://github.com/ friendly/matlib/issues/58#issuecomment-2303586612, or unsubscribe https://github.com/notifications/unsubscribe-auth/ ADLSAQU67XI3LYXXCNX4KBTZSVIWFAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBTGU4DMNRRGI. You are receiving this because you are subscribed to this thread.Message ID: @.***>

— Reply to this email directly, view it on GitHubhttps://github.com/friendly/matlib/issues/58#issuecomment-2304434726, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AADPNIKKQEG6KQSQ3NO2ZE3ZSXDN3AVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBUGQZTINZSGY. You are receiving this because you authored the thread.Message ID: @.***>

john-d-fox commented 2 months ago

On 2024-08-22 7:38 a.m., Michael Friendly wrote:

Caution: External email.

Another option is to have a parallel function, matmult(A, B), that takes a simplify argument More on this later

I don't see a need for that, and simplification is more general than for matrix multiplication.

Because I can't get back to this until later in the day, I've uploaded simplify.R to dev/ with the work I've done so far -- a generic simplify() function with a "default" method and some examples. I need to add a "latexMatrix" method and do some more testing.

John

Get Outlook for iOShttps://aka.ms/o0ukef


From: John Fox @.> Sent: Thursday, August 22, 2024 1:28:29 PM To: friendly/matlib @.> Cc: Michael L Friendly @.>; Author @.> Subject: Re: [friendly/matlib] latexMatrix arithmetic: trivial simplifications? (Issue #58)

I think that sparse = TRUE is different: replacing 0s with blanks. That's a worthwhile option, but it's different from, and simpler than, replacing expressions that evaluate to 0 with 0 (which the kronecker() method also does) and removing multiplication by 1.

I made good progress on a simplify() function yesterday evening. I have appointments today, so probably won't get much done, but should have a preliminary implementation soon.

simplify() could optionally be called by print.latexMatrix() or used directly.

John

On 2024-08-21 11:07 p.m., Phil Chalmers wrote:

Caution: External email.

It's certainly possible to do this. I did something similar in |dev/ Arith2.R| for the kronecker product, and I was glad to see this idea made into John's implementation as well. Perhaps the best way to do this is to add a |print(..., sparse = TRUE)| to detect any zero multiplications in defined matrices and have them drop out?

— Reply to this email directly, view it on GitHub https://github.com/ friendly/matlib/issues/58#issuecomment-2303586612, or unsubscribe <https://github.com/notifications/unsubscribe-auth/

ADLSAQU67XI3LYXXCNX4KBTZSVIWFAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBTGU4DMNRRGI>. You are receiving this because you are subscribed to this thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub<https://github.com/ friendly/matlib/issues/58#issuecomment-2304434726>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ AADPNIKKQEG6KQSQ3NO2ZE3ZSXDN3AVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBUGQZTINZSGY>. You are receiving this because you authored the thread.Message ID: @.***>

— Reply to this email directly, view it on GitHub <https://github.com/ friendly/matlib/issues/58#issuecomment-2304452996>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ ADLSAQRLJTNFGBTSNQRDCOTZSXETRAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBUGQ2TEOJZGY>. You are receiving this because you commented.Message ID: <friendly/ @.***>

john-d-fox commented 2 months ago

I had some time between appointments, so I added a simplify.latexMatrix() method to dev/simplify.R.

For Michael's original example, this produces

> simplify(C %*% B)
\begin{bmatrix}  
\beta_{1,y_{1}} & \beta_{1,y_{2}} & \beta_{1,y_{3}} \\ 
\beta_{2,y_{1}} & \beta_{2,y_{2}} & \beta_{2,y_{3}} \\ 
\end{bmatrix}

In my opinion, that's not more complicated than something like matmult(C, B, simplify=TRUE).

The code in dev/simplify.R could use more testing, and I expect that it could also be improved -- I'm not a regular expression whiz.

friendly commented 2 months ago

simplify.R is very nice. Hard to believe it was something you could do between appointments!

But I was thinking of some way, hopefully more elegant, that would avoid having to parse what latexMatrix generates. That is, the best place to simplify would be as source, where the elements of the dot product, dot(a, b) can be recognized as in 3 cases:

Just for fun, I'll try to sketch out a function that works this way.

john-d-fox commented 2 months ago

I added a matmult() function to dev/simplify.R. The rationale is that the function takes an arbitrary number of "latexMatrix" objects as arguments and multiplies them in left-to-right order. There are examples in the file. All of the following work, e.g.:

matmult(R)
matmult(R, S)
matmult(R, S, T)

assuming that R, S, and T are conformable for multiplication in the order given.

john-d-fox commented 2 months ago

Hi,

On 2024-08-22 2:43 p.m., Michael Friendly wrote:

Caution: External email.

|simplify.R| is very nice. Hard to believe it was something you could do between appointments!

I just wrote the "latexMatrix" method between appointments; that's pretty simple.

But I was thinking of some way, hopefully more elegant, that would avoid having to parse what |latexMatrix| generates. That is, the best place to simplify would be as source, where the elements of the dot product, |dot(a, b)| can be recognized as in 3 cases:

  • 0 in either a[i ] or b[j-]> ignore term
  • 1 -> use the other
  • -1 -> use the other, and combine with |-| For another thing, |\cdot| is hard-coded throughout, and we might at sometime want to change it or allow an argument, `times = c("\cdot", "\times, "")

Just for fun, I'll try to sketch out a function that works this way.

— Reply to this email directly, view it on GitHub <https://github.com/ friendly/matlib/issues/58#issuecomment-2305411510>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ ADLSAQQVVN52ZWXZLRYZF53ZSYWN5AVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBVGQYTCNJRGA>. You are receiving this because you commented.Message ID: <friendly/ @.***>

friendly commented 2 months ago

I added dev/dot.R, doing the latex dot product of two vectors, either numeric or symbolic, and simplifying 0s & 1s. This avoids the need to parse the results as in simplify()

>  num <- 0:2
>  chr <- letters[1:3]
>  
>  dot(num, num)
[1] "2 \\cdot 2"
>  dot(num, chr)
[1] "b + 2 \\cdot c"
>  dot(chr, num)
[1] "b + c \\cdot 2"
>  dot(chr, chr)
[1] "a \\cdot a + b \\cdot b + c \\cdot c"

If it was OK to just simplify %*% by default (I see nothing wrong with that), this could be used inside the loop at lines 724 - 730 of latexMatrix.R

      for (k in 1:ncol(X)){
        Z[i, j] <- paste0(Z[i, j],
                          if (k > 1) " + ",
                          parenthesize(X[i, k]),
                          " \\cdot ",
                          parenthesize(Y[k, j]))
      }

replacing this with

Z[i,j] <- dot(X[i,], y[,j]
john-d-fox commented 2 months ago

I tried this, fixing the typos in the last line: Z[i, j] <- dot(X[i, ], Y[, j])

The resulting %*% fails for the following example (the second one that I tested):

> B <- latexMatrix('\\beta', ncol = 3, nrow=4, 
+                  comma=TRUE, prefix.col = 'y_',
+                  zero.based=c(TRUE, FALSE))
> D <- diag(4)
> D[3, 3] <- 0
> D <- latexMatrix(D)
> D %*% B
\begin{pmatrix}  
\beta_{0,y_{1}} & \beta_{0,y_{2}} & \beta_{0,y_{3}} \\ 
\beta_{1,y_{1}} & \beta_{1,y_{2}} & \beta_{1,y_{3}} \\ 
                &                 &                 \\ 
\beta_{3,y_{1}} & \beta_{3,y_{2}} & \beta_{3,y_{3}} \\ 
\end{pmatrix}

Some comments:

I like your approach, which I'm sure can be made more robust. As you suggest, it's simpler than simplifying at the end. You need to test with a more diverse set of examples and deal with cases like the one above.

The reason that I originally didn't simplify by default is that I thought that it was useful to show how the matrix product is formed, given the general purpose of the matlib package. That could still be done by including the matmult() function, with optional simplification, along with the %*% operator.

You mentioned that you'd like to make the LaTeX multiplication character flexible. I think that can be done by supplying a variable in the package, say latexTimesChar set initially to "\\cdot", along with a mechanism for changing it.

john-d-fox commented 2 months ago

On second thought, a latexTimesChar variable is more complicated than an option(), and still violates functional programming. I've gone ahead and implemented a latexMultSymbol option(); if it's not set, "\\cdot" is used.

friendly commented 2 months ago

I see why your %*% example failed, where the 0 is actually desired. This calls for a simplify option, but not sure where to implement it. That's the trouble with operators, where you can't have optional arguments. Except, a hint:

get('%*%')(A, B, simplify=TRUE)

Or in matmult()

john-d-fox commented 2 months ago

I think that makes the issue more complicated than it needs to be. You should be able to handle this and other problems that may arise in dot() with more extensive testing.

In this case, adding the line

  if (res == "") res <- "0"

as the next-to-last command in dot() produces

> D %*% B
\begin{pmatrix}  
\beta_{0,y_{1}} & \beta_{0,y_{2}} & \beta_{0,y_{3}} \\ 
\beta_{1,y_{1}} & \beta_{1,y_{2}} & \beta_{1,y_{3}} \\ 
0               & 0               & 0               \\ 
\beta_{3,y_{1}} & \beta_{3,y_{2}} & \beta_{3,y_{3}} \\ 
\end{pmatrix}

as desired.

More generally, however, I think that there's an argument for making matmult() the place where the computations happen, which would allow for additional arguments like simplify, and then to have %*% call matmult() (rather than vice-versa) with arguments set to specific values like simplify = TRUE.

Since I spent quite a bit of time trying to make simplify() robust, I'd rather not do the same for dot(), but I expect that the latter will end up being considerably simpler than the former. That is, I think that your approach is better, though it now requires additional work to make it robust.

john-d-fox commented 2 months ago

Another issue -- handling -1s:

> F <- latexMatrix(matrix(c(-1, 0, 1, 0, 0, -1, 0, 1), ncol=4))
> F
\begin{pmatrix} 
-1 &  1 &  0 &  0 \\ 
 0 &  0 & -1 &  1 \\ 
\end{pmatrix}
> F %*% B
\begin{pmatrix}  
(-1) \cdot \beta_{0,y_{1}} + \beta_{1,y_{1}} & (-1) \cdot \beta_{0,y_{2}} + \beta_{1,y_{2}} & (-1) \cdot \beta_{0,y_{3}} + \beta_{1,y_{3}} \\ 
(-1) \cdot \beta_{2,y_{1}} + \beta_{3,y_{1}} & (-1) \cdot \beta_{2,y_{2}} + \beta_{3,y_{2}} & (-1) \cdot \beta_{2,y_{3}} + \beta_{3,y_{3}} \\ 
\end{pmatrix}
friendly commented 2 months ago

Yes, I had that on my todo list, but didn't do it in my initial draft, partly because that was a little trickier given the initial logic.

friendly commented 2 months ago

I generally agree with your earlier comments about simplify and matmult, and I think you're in a better position to decide how to handle something like dot in this mix. I'll add your '0' fix and perhaps we can both identify where/how it can be made more robust.

When would it make sense to move simplify.R to R/?

john-d-fox commented 2 months ago

I was about to propose something very similar, so I think we're on the same wavelength.

I suggest the following:

  1. You work to make dot() more robust. At this point, the only issues I've uncovered are replacing blanks with 0s, which as you note is easily fixed, and handling multiplication by -1. You can add a simplify argument, which defaults to TRUE to dot().
  2. When you're done, I rewrite matmult() so that it uses dot() to perform matrix multiplication. It too will have a simplify argument defaulting to TRUE, which is just passed to dot().
  3. I rewrite the %*%.matrixMult() method so that it calls matmult().
  4. I take another look at the other arithmetic functions and operators to see whether a similar approach to simplification is desirable and feasible.

I don't think that it's a good idea to move simplify.R to R/. The code for simplify() is kludgy and convoluted (which is why I prefer your dot()-based approach), and it should be unnecessary once the above changes are implemented.

Does all that sound reasonable?

john-d-fox commented 2 months ago

I hope that this isn't counterproductive, but for curiosity, I tried to fix dot() and came up with the following:

dot <- function(x, y) {
  if (length(x) != length(y)) stop("Vectors must have the same length")
  x <- trimws(x)
  y <- trimws(y)
  res <- ""
  for (i in 1:length(x)) {
    # ignore terms multiplied by zero
    if (x[i] == "0") next
    if (y[i] == "0") next
    xi <- if(x[i] == "1") "" else x[i]
    yi <- if(y[i] == "1") "" else y[i]
    if (x[i] == "1" && y[i] == "1") xi <- "1"
    xi <- if (xi == "-1") "-" else xi
    if (y[i] == "-1") {
      yi <- ""
      xi <- if (x[i] == "-1") "1" else paste0("-", parenthesize(xi))
    }
    times <- if(xi == "" || xi == "-" || yi == "") "" else " \\cdot "
    res <- paste0(res,
                  if (nchar(res) > 0) " + ",
                  if (y[i] != "-1" && xi != "-") parenthesize(xi) else xi,
                  times,
                  parenthesize(yi))
  }
  if (res == "") res <- "0"
  res
}

it seems to work OK.

friendly commented 2 months ago

OK, I'll work from your version of dot. It should take one simplify=T/F arg I think not separate ones for 1s & 0s.

In terms of making it 'robust', I'm not sure what to test for, other than the combinations of num and chr, now including -1, and also with simplify = FALSE

friendly commented 2 months ago

OK, I incorporated simplify = T/F in a new dev/dot.R. Seems to work for the tests I've run

>  # w/ simplify=T
>   dot(num, num)  # OK
[1] "1 + 1 + 2 \\times 2"
>   dot(num, chr)  # OK
[1] "-a + c + 2 \\times d"
>   dot(chr, num)  # OK
[1] "-a + c + d \\times 2"
>   dot(chr, chr)  # OK
[1] "a \\times a + b \\times b + c \\times c + d \\times d"

W/o simplify

>   dot(num, num, simplify = FALSE) 
[1] "(-1) \\times (-1) + 0 \\times 0 + 1 \\times 1 + 2 \\times 2"
>   dot(num, chr, simplify = FALSE)
[1] "(-1) \\times a + 0 \\times b + 1 \\times c + 2 \\times d"
>   dot(chr, num, simplify = FALSE)
[1] "a \\times (-1) + b \\times 0 + c \\times 1 + d \\times 2"
>   dot(chr, chr, simplify = FALSE)
[1] "a \\times a + b \\times b + c \\times c + d \\times d"
>   # change the mult symbol  
>   opt <- options(latexMultSymbol = "\\times") 
>   options("latexMultSymbol")
$latexMultSymbol
[1] "\\times"

>   dot(num, num)
[1] "1 + 1 + 2 \\times 2"
>   dot(num, chr)
[1] "-a + c + 2 \\times d"
friendly commented 2 months ago

Sorry; I copy/pasted those tests out of order, after I had changed latexMultSymbol to \\times

>   options(latexMultSymbol = "\\cdot")
>   
>   dot(num, num)  # OK
[1] "1 + 1 + 2 \\cdot 2"
>   dot(num, chr)  # OK
[1] "-a + c + 2 \\cdot d"
>   dot(chr, num)  # OK
[1] "-a + c + d \\cdot 2"
>   dot(chr, chr)  # OK
[1] "a \\cdot a + b \\cdot b + c \\cdot c + d \\cdot d"
john-d-fox commented 2 months ago

All this looks promising. With respect to additional testing, you can create some matrices and test with %*% using more complex matrix expressions like A %*% (B - C), A %*% -B, (A - k*D) %*% B (where k is a scalar), A %*% B %*% E, A %*% solve(B), etc., and see whether the simplifications hold up.

friendly commented 2 months ago

Then, I'd have to incorporate dot() into latexMartrix(). I suppose I can do this from a copy of the latter in dev/.

john-d-fox commented 2 months ago

Here's what I'm using:

numericDimensions <- matlib:::numericDimensions
updateWrapper <- matlib:::updateWrapper
parenthesize <- matlib:::parenthesize

`%*%.latexMatrix` <- function(x, y){
  if (!inherits(y, "latexMatrix")){
    stop(deparse(substitute(y)),
         " is not of class 'latexMatrix'")
  }
  numericDimensions(x)
  numericDimensions(y)
  X <- getBody(x)
  Y <- getBody(y)
  dimX <- dim(X)
  dimY <- dim(Y)
  if (dimX[2] != dimY[1]){
    stop('matricies are not conformable for multiplication')
  }

  Z <- matrix("", nrow(X), ncol(Y))
  for (i in 1:nrow(X)){
    for (j in 1:ncol(Y)){
      for (k in 1:ncol(X)){
        Z[i, j] <- dot(X[i, ], Y[, j])
      }
    }
  }
  result <- latexMatrix(Z)
  result <- updateWrapper(result, getWrapper(x))
  result$dim <- dim(Z)
  result
}
friendly commented 2 months ago

OK -- your test code above, along with some tests is now in dev/dot-test.R.

My first test gives a strange result, which I don't see the reason for

 (A <- latexMatrix(matrix(c(1, -3, 0, 1), 2, 2)))
 (B <- latexMatrix(matrix(c(5, 3, -1, 4), 2, 2)))
 A %*% B

Gives a naked - sign n row 1, col 2.

\begin{pmatrix}  
5                & -         \\ 
(-3) \cdot 5 + 3 & -(-3) + 4 \\ 
\end{pmatrix}

We're going out to dinner shortly, so I'm probably done for the day.

john-d-fox commented 2 months ago

A simple fix is to add if (res == "-") res <- "-1" near the bottom of dot(), resulting in

\begin{pmatrix}  
5                & -1        \\ 
(-3) \cdot 5 + 3 & -(-3) + 4 \\ 
\end{pmatrix}

That's correct but clearly could be simplified further. BTW, for strictly numeric matrices, one could always do the numeric matrix multiplication and then make the "latexMatrix" object.

john-d-fox commented 2 months ago

Here's another approach, which converts the matrices to numeric if that's possible.

Code:

`%*%.latexMatrix` <- function(x, y){
  if (!inherits(y, "latexMatrix")){
    stop(deparse(substitute(y)),
         " is not of class 'latexMatrix'")
  }
  numericDimensions(x)
  numericDimensions(y)
  X <- getBody(x)
  Y <- getBody(y)
  dimX <- dim(X)
  dimY <- dim(Y)
  if (dimX[2] != dimY[1]){
    stop('matricies are not conformable for multiplication')
  }

  XX <- suppressWarnings(as.numeric(X))
  YY <- suppressWarnings(as.numeric(Y))
  if (!any(is.na(XX)) && !any(is.na(YY))){
    XX <- matrix(XX, dimX[1], dimX[2])
    YY <- matrix(YY, dimY[1], dimY[2])
    Z <- XX %*% YY

  } else {

    Z <- matrix("", nrow(X), ncol(Y))
    for (i in 1:nrow(X)){
      for (j in 1:ncol(Y)){
        for (k in 1:ncol(X)){
          Z[i, j] <- dot(X[i, ], Y[, j])
        }
      }
    }
  }

  result <- latexMatrix(Z)
  result <- updateWrapper(result, getWrapper(x))
  result$dim <- dim(Z)
  result
}

dot <- function(x, y) {
  if (length(x) != length(y)) stop("Vectors must have the same length")
  x <- trimws(x)
  y <- trimws(y)
  res <- ""
  for (i in 1:length(x)) {
    # ignore terms multiplied by zero
    if (x[i] == "0") next
    if (y[i] == "0") next
    xi <- if(x[i] == "1") "" else x[i]
    yi <- if(y[i] == "1") "" else y[i]
    if (x[i] == "1" && y[i] == "1") xi <- "1"
    xi <- if (xi == "-1") "-" else xi
    if (y[i] == "-1") {
      yi <- ""
      xi <- if (x[i] == "-1") "1" else paste0("-", parenthesize(xi))
    }
    times <- if(xi == "" || xi == "-" || yi == "") "" else " \\cdot "
    res <- paste0(res,
                  if (nchar(res) > 0) " + ",
                  if (y[i] != "-1" && xi != "-") parenthesize(xi) else xi,
                  times,
                  parenthesize(yi))
  }
  if (res == "") res <- "0"
  if (res == "-") res <- "-1"
  res
}

Behaviour:

> C <- latexMatrix(matrix(c(0, 1, 0, 0,
+                           0, 0, 1, 0), nrow=2, byrow=TRUE), 
+                  matrix = "bmatrix")
> B <- latexMatrix('\\beta', ncol = 3, nrow=4, 
+                  comma=TRUE, prefix.col = 'y_',
+                  zero.based=c(TRUE, FALSE))
> C %*% B
\begin{bmatrix}  
\beta_{1,y_{1}} & \beta_{1,y_{2}} & \beta_{1,y_{3}} \\ 
\beta_{2,y_{1}} & \beta_{2,y_{2}} & \beta_{2,y_{3}} \\ 
\end{bmatrix}

> A <- latexMatrix(matrix(c(1, -3, 0, 1), 2, 2))
> B <- latexMatrix(matrix(c(5, 3, -1, 4), 2, 2))
> A %*% B
\begin{pmatrix}  
  5 &  -1 \\ 
-12 &   7 \\ 
\end{pmatrix}
friendly commented 2 months ago

I like this very much. I don't see any reason to stick with latex expressions that evaluate to a number when both X and Y are numeric. Needs more testing, though

friendly commented 2 months ago

Your code and a bunch of tests is now in dev/dot-test2.R.

I'm going to move dev/dot.R to R/ and finish documenting it. For now in a separate .Rd. The docs for latexMatrix are now quite unwieldly. I think that operation / operators for latex matrices would be better in a separate .Rd file.

john-d-fox commented 2 months ago

I agree that it would be best to document the arithmetic functions and operators separately. Also, defaulting to numeric operations on numeric matrices can be done in all of the other matrix arithmetic functions and operators. I think that it would be useful to have arithmetic functions corresponding to the other operators, like matmult() for %*% -- add() for +, etc. -- with the actual computations done in the functions, which each would have a simplify operator defaulting to TRUE. The operators would then just call the corresponding functions with default arguments. That would allow, e.g., showing details when desired but suppressing them by default. (I made a similar suggestions previously.) I likely won't have much time today to work on this but should have plenty of time during the week.

friendly commented 2 months ago

We agree on your first two statements. I'm not entirely clear what you have in mind for the other matrix operations.

matmult.latexMatrix <- function(X, ..., simplify=TRUE){
  matrices <- list(...)
  for (i in seq_along(matrices)){
    X <- X %*% matrices[[i]]
  }
}
john-d-fox commented 2 months ago

Sorry I wasn't clear. The basic idea is that the computations would, e.g. be done in matmult(), which, in addition to being able to handle more than two matrices, would have simplify and as.numeric arguments, which default to TRUE, but which could be FALSE, thus showing the details of the cell products. (The as.numeric argument would if TRUE convert the matrices to numeric to get maximum simplification prior to returning a "latexMatrix" object, as in the latest version of %*%.) %*% would then just call matmult(), rather than vice-versa, with two matrix arguments and (implicitly) simplify=TRUE and as.numeric=TRUE. Similarly, add(), not +, would do the actual computation. Like matmult(), it could handle more than two matrices and have simplify and as.numeric arguments defaulting to TRUE. + would call add() with two matrix arguments and simplify=TRUE and as.numeric=TRUE. Likewise for other matrix operators. As I said, I made this suggestion previously. It should be easy to implement and has the advantage of maximum flexibility, provided by the mathematical functions, and convenient default behaviour, provided by the operators.

friendly commented 2 months ago

OK, that is clear, and it sounds like a good plan. I think it makes the most sense for me not to change anything until you have a chance to sketch this out in code.

It might make sense to put these in a file, matops.R or matrixOps.R along with my dot.R, which could be documented together with the title "Operations on LaTeX Matrices` or something like that.

john-d-fox commented 2 months ago

On 2024-08-25 11:43 a.m., Michael Friendly wrote:

Caution: External email.

OK, that is clear, and it sounds like a good plan. I think it makes the most sense for me not to change anything until you have a chance to sketch this out in code.

That makes sense.

It might make sense to put these in a file, |matops.R| or |matrixOps.R| along with my |dot.R|, which could be documented together with the title "Operations on LaTeX Matrices` or something like that.

That's pretty much what I had in mind.

— Reply to this email directly, view it on GitHub https://github.com/friendly/matlib/issues/58#issuecomment-2308901042, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADLSAQVW7NTHGKEYEX4IYXDZTH3RDAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMBYHEYDCMBUGI. You are receiving this because you commented.Message ID: @.***>

friendly commented 2 months ago

I was looking at kronecker(), still in dev/. I wondered if there was a way, with latexMatrix and friends to generate the "definition" of Kronecker product: each element of A & all of B:

\mathbf{A} \otimes \mathbf{B} = &
\begin{pmatrix}
a_{11} \cdot \mathbf{B} & a_{12} \cdot \mathbf{B} \\
a_{21} \cdot \mathbf{B} & a_{22} \cdot \mathbf{B}
\end{pmatrix} 

image

john-d-fox commented 2 months ago

The most current code for the kronecker() method is in R/latexMatrix.R, at least until I move it to latexMatrixOperations.R. It may be the same as the earlier code in dev/kronecker.R -- I haven't compared the two. As to the question: I'm sure that one could do what you want but the behaviour would be different from every other "latexMatrix" function and operator. If you want to give this a try, introducing say an argument expand=FALSE, you can go ahead. Alternatively, we might want to support partitioned matrices more directly (if simply bolding matrix elements doesn't already work -- not tried). At the moment, I'm busy with writing some new functions, like matsum() and matpower, and dividing functions and operators between latexMatrix.R and latexMatrixOperations.R.

philchalmers commented 2 months ago

Hi Michael,

Wouldn't the definition just be this?

A <- latexMatrix('a', nrow=2, ncol=2)
B <- latexMatrix('\\mathbf{B}', ncol=1, nrow=1)
kronecker(A, B)

\begin{pmatrix}  
a_{11} \cdot \mathbf{B} & a_{12} \cdot \mathbf{B} \\ 
a_{21} \cdot \mathbf{B} & a_{22} \cdot \mathbf{B} \\ 
\end{pmatrix}
friendly commented 2 months ago

@john-d-fox My Q about kronecker() was perhaps out of place in this thread. Phi's reply was sufficient for what I was seeking. I didn't mean to divert you or suggest that we modify it to handle partitioned matrices.

john-d-fox commented 2 months ago

Phil's clever example shows how to handle partitioned matrices without any modifications to the current code! (I thought something like that might work, which is why I said "if simply bolding matrix elements doesn't already work -- not tried.")

On 2024-08-27 1:54 p.m., Michael Friendly wrote:

Caution: External email.

@john-d-fox https://github.com/john-d-fox My Q about |kronecker()| was perhaps out of place in this thread. Phi's reply was sufficient for what I was seeking. I didn't mean to divert you or suggest that we modify it to handle partitioned matrices.

— Reply to this email directly, view it on GitHub <https://github.com/ friendly/matlib/issues/58#issuecomment-2313183222>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ ADLSAQSOPX65HYZ67FK3THDZTS4MLAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMJTGE4DGMRSGI>. You are receiving this because you were mentioned.Message ID: <friendly/ @.***>

john-d-fox commented 2 months ago

I just pushed a new version of R/latexMatrix.R along with the new file R/latexMatrixOperations.R. These require some more work on the Roxygen markup, but the package builds and checks without errors or notes.

Two problems that I've not been able to solve:

(1) I wanted the .Rd file for the arithmetic functions and operators to be latexMatrixOperations.Rd, but I don't know how to do that in Roxygen since latexMatrixOperations isn't a function name. At present, the file name is matsum.Rd, the first function in the file. I haven't tried very hard to figure out how to do this, thinking one of you might know.

(2) I've marked a couple of examples in latexMatrixOperations.R as \dontrun{} because not doing that caused an error in R CMD check even though the examples run perfectly well if they're included in the .Rd file.

friendly commented 2 months ago

Just to complete Phil's example (perhaps for the vignette):

  A <- latexMatrix("a", 2, 2)
  B <- latexMatrix("b", 2, 2)
  kronecker(A, B)

  # Generate the 'definition' of Kronecker product,
  Bmat <- latexMatrix('\\mathbf{B}', ncol=1, nrow=1)
  kronecker(A, Bmat)

  Eqn("\\mathbf{A} \\otimes \\mathbf{B} = &",
      kronecker(A, Bmat),
      "\\\\[1.5ex]\n= & ",
      kronecker(A, B),
      align = TRUE)

image

john-d-fox commented 2 months ago

OK, I quickly figured out how to use latexMatrixOperations.Rd (problem 1 identified above).

john-d-fox commented 2 months ago

Problem 2 (examples marked as \dontrun{}) now fixed. Source of the problem was dot.R changing the LaTeX multiplication character and \cdothard-coded in as.double.latexMatrix(). dot.R now moved to dev/ (it was redundant) and \cdot no longer hard-coded, including in showEqn.R.

friendly commented 2 months ago

Can’t wait to see. My internet 🛜 is down

Get Outlook for iOShttps://aka.ms/o0ukef


From: John Fox @.> Sent: Tuesday, August 27, 2024 10:58:41 PM To: friendly/matlib @.> Cc: Michael L Friendly @.>; Author @.> Subject: Re: [friendly/matlib] latexMatrix arithmetic: trivial simplifications? (Issue #58)

Problem 2 (examples marked as \dontrun{}) now fixed. Source of the problem was dot.R changing the LaTeX multiplication character and \cdot hard-coded in as.double.latexMatrix(). dot.R now moved to dev/ (it was redundant) and \cdot no longer hard-coded, including in showEqn.R.

— Reply to this email directly, view it on GitHubhttps://github.com/friendly/matlib/issues/58#issuecomment-2313526064, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AADPNIMTQWMG63KGOOBFUJTZTTSADAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMJTGUZDMMBWGQ. You are receiving this because you authored the thread.Message ID: @.***>

friendly commented 2 months ago

I edited the documentation of latexMatrixOperations a bit. This looks really nice.

I'm beginning to extend the latex-equations vignette, and would like to do something on partitioned matrices, e.g.,a 4×4 matrix $M$, which is partitioned in 2×2 b locks, which are labeled $M{i,j}$.

M <- latexMatrix("m", 4, 4)
Mpart <- latexMatrix('\\mathbf{M}', nrow = 2, ncol = 2, comma = TRUE)
Eqn("\\mathbf{M} =", Mpart,
    " =", M)

Perhaps it would be useful to implement [.latexMatrix() to select rows / cols?

john-d-fox commented 2 months ago

Please see the issue "conflicts in latexMatrixOperations.R," which I filed about an hour ago.

I can implement [.latexMatrix() as long as you're not going to work on the file latexMatrix.R (to avoid additional conflicts). As I said in the issue I filed, I'll now stay away from latexMatrixOperations.R.

John

On 2024-08-28 1:10 p.m., Michael Friendly wrote:

Caution: External email.

I edited the documentation of |latexMatrixOperations| a bit. This looks really nice.

I'm beginning to extend the |latex-equations| vignette, and would like to do something on partitioned matrices, e.g.,a 4×4 matrix $M$, which is partitioned in 2×2 b locks, which are labeled $M{i,j}$.

|M <- latexMatrix("m", 4, 4) Mpart <- latexMatrix('\mathbf{M}', nrow = 2, ncol = 2, comma = TRUE) Eqn("\mathbf{M} =", Mpart, " =", M) |

Perhaps it would be useful to implement |[.latexMatrix()| to select rows / cols?

— Reply to this email directly, view it on GitHub <https://github.com/ friendly/matlib/issues/58#issuecomment-2315869972>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ ADLSAQTEQYORR5IELTQZCJDZTX77PAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMJVHA3DSOJXGI>. You are receiving this because you were mentioned.Message ID: <friendly/ @.***>

john-d-fox commented 2 months ago

Indexing turns out to be simple:

`[.latexMatrix` <- function(x, i, j, ..., drop){
  matlib:::numericDimensions(x)
  X <- getBody(x)
  X <- X[i, j, drop=FALSE]
  X <- latexMatrix(X)
  matlib:::updateWrapper(X, getWrapper(x))
}

Of course, matlib::: references can be removed once this is in the package. (I'm waiting to do that pending confirmation that no one is currently working on latexMatrix.R.)

philchalmers commented 2 months ago

Beat me to it :) I also had plans to add rbind() and cbind() support, mainly for convenient presentation purposes (e.g., building partitioned matrices, like these https://tex.stackexchange.com/questions/262426/how-to-create-a-simple-partitioned-matrix). Though I likely won't have time to work on these until elementary school restarts...

john-d-fox commented 2 months ago

That's also straightforward. The following seems to work:

cbind.latexMatrix <- function(..., deparse.level = 1){
  matrices <- list(...)
  nrow <- Nrow(matrices[[1]])
  if (length(matrices) == 1) return(matrices[[1]])
  else {
    X <- getBody(matrices[[1]])
    for (matrix in matrices[-1]){
      if (Nrow(matrix) != nrow) stop("number of rows not the same")
      Y <- getBody(matrix)
      X <- cbind(X, Y)
    }
    X <- latexMatrix(X)
    return(matlib:::updateWrapper(X, getWrapper(matrices[[1]])))
  }
}

and similarly for rbind().

I've put [.latexMatrix(), cbind.latexMatrix(), and rbind.latexMatrix() in dev/bracket.R, along with some examples.

If this looks OK, and no one else is modifying R/latexMatrix.R, I can move these functions there.

friendly commented 2 months ago

So sorry for the mixup. I thought you were done w/ latexMatrixOps for now. I have an appt tomorrow so won’t be doing any work til late in my day

Get Outlook for iOShttps://aka.ms/o0ukef


From: John Fox @.> Sent: Wednesday, August 28, 2024 10:11:05 PM To: friendly/matlib @.> Cc: Michael L Friendly @.>; Author @.> Subject: Re: [friendly/matlib] latexMatrix arithmetic: trivial simplifications? (Issue #58)

That's also straightforward. The following seems to work:

cbind.latexMatrix <- function(..., deparse.level = 1){ matrices <- list(...) nrow <- Nrow(matrices[[1]]) if (length(matrices) == 1) return(matrices[[1]]) else { X <- getBody(matrices[[1]]) for (matrix in matrices[-1]){ if (Nrow(matrix) != nrow) stop("number of rows not the same") Y <- getBody(matrix) X <- cbind(X, Y) } X <- latexMatrix(X) return(matlib:::updateWrapper(X, getWrapper(matrices[[1]]))) } }

and similarly for rbind().

I've put [.latexMatrix(), cbind.latexMatrix(), and rbind.latexMatrix() in dev/bracket.R, along with some examples.

If this looks OK, and no one else is modifying R/latexMatrix.R, I can move these functions there.

— Reply to this email directly, view it on GitHubhttps://github.com/friendly/matlib/issues/58#issuecomment-2316164821, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AADPNII7Z2G6AW43DAIDFYTZTYVFTAVCNFSM6AAAAABM42L7RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMJWGE3DIOBSGE. You are receiving this because you authored the thread.Message ID: @.***>