Msg from Brian Ripley:

R-devel has a new runtime check in sprintf() that all arguments are used by the format and this is showing one or more warnings in the checks of these packages, usually in the examples. E.g.

Warning in sprintf("%s$", extension, ignore.case = TRUE) :
argument not used by format
Warning in sprintf("%s", getFunctionLinks(tokens[w]), tokens[w], :
argument not used by format

Please investigate (by running R CMD check with the latest R-devel and grepping the .Rcheck directory). In some cases the warning seems to be from another package, usually one with the same maintainer. If from one not on this list, please let CRAN and its maintainer know. If you think this is a false positive, please contact Tomas Kalibera (in Cc:).

At present, there is no official NEWS file in the package appearing on CRAN or available via news("matlib"). The devtools / Github way seems to be to just have a file NEWS.md in the
top level of the package, but this is not recognized in CRAN builds.

In the Makefile, I now have a target "news", so that make news does a minimal conversion to
a NEWS file. Is there a better way?

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

Ax = b
LUx = b
LUx - b = 0
L(Ux - d) = 0

where Ld = b. Then do a simple forward-solve for d, and finally a simple back-solve for x given d.

Any interest in this form of decomposition in the form of a function LU()? With all the ERO already defined it should be easy to do, and the solving of d/x could be made into a flag solve = TRUE. Otherwise, just the L and U matrices will be returned.

In this query to stackoverflow, http://stackoverflow.com/questions/39419622/printing-matrices-and-vectors-side-by-side, I ask about code to print a collection of matrices, vectors, and character strings side-by-side, to illustrate matrix equations and expressions. My main example of what I want to do is this expression from SAS/IML

print A ' * ' x '=' (A * x) '=' b;

                A                     X   #TEM1001       B
                1     1    -4  *  0.733 =        2 =     2
                1    -2     1     -0.33          1       1
                1     1     1      -0.4          0       0

There was one answer so far, but not as general as I'd like. Is it posible to do better?

Hej,
the matlib::echelon function cannot handle matrices with one row. While this may be usually not needed at all, I am using the echelon function in a package where it is possible that sometimes the A matrix can indeed have one row only. I would prefer the function not to break the code run with an error.

A <- matrix(c(0,1,1,0), 2)
b <- 1:2
matlib::echelon(A,b) # works
matlib::echelon(A[1, , drop = F],b[1]) #fails
# should return sth like
#      [,1] [,2]  [,3]                                                                                    
  [1,]    0    1    1

John & I discussed some issues relating to development of the package toward a next release:

• Naming conventions: At present, we use a mix of camelCase (gaussianElimination, showEqn, ...) and underscored_names (power_method). We should move toward a common standard, e.g., power_method -> powerMethod. This should probably apply to file names also: e.g., gaussianElimination() is in gaussian-elimination.R.
• There are a few functions that are redundant (e.g, echelon is a call to gaussianElimination) or synonyms (inv is a call to Inverse), but where they serve a semantic purpose, we should allow this limited use.
• Where we supersede a standard function in R, we should use a consistent convention, e.g., an initial upper case letter, unless the function name is an acronym. Thus, qr() -> QR(), and svd() -> SVD() are OK, but eigen() should be renamed Eigen(), rather than eig().
• This points to the usefulness of a Det() function that would extend det() with a verbose argument.
• We should use our own functions consistently in other functions.

Hi,

I thought about adding a GramSchmidt() function and discovered that there already is one! The problem is that it can easily fail misleadingly if the matrix isn't of full column rank. Here's an improved version, and an example -- try the example first with the current version to see what I mean:

GramSchmidt <- function (X, normalize = TRUE, verbose = FALSE, 
                         tol=sqrt(.Machine$double.eps)) {
    B <- X
    B[, 2L:ncol(B)] <- 0
    if (verbose) {
        cat("\nInitial matrix:\n")
        printMatrix(X)
        cat("\nColumn 1: z1 = x1\n")
        printMatrix(B)
    }
    for (i in 2:ncol(X)) {
        res <- X[, i]
        for (j in 1:(i - 1)) res <- res - Proj(X[, i], B[, j])
        if (len(res) < tol) res <- 0
        B[, i] <- res
        if (verbose) {
            prj_char <- character(i - 1)
            for (j in 1:(i - 1)) prj_char[j] <- sprintf("Proj(x%i, z%i)", 
                                                        i, j)
            cat(sprintf("\nColumn %i: z%i = x%i - %s\n", i, i, 
                        i, paste0(prj_char, collapse = " - ")))
            printMatrix(B)
        }
    }
    B <- B[, !apply(B, 2, function(b) all(b == 0))]
    if (normalize) {
        norm <- diag(1/len(B))
        B <- B %*% norm
        if (verbose) {
            cat("\nNormalized matrix: Z * inv(L) \n")
            printMatrix(B)
        }
    }
    if (verbose) 
        return(invisible(B))
    B
}

set.seed(123)
X <- matrix(rnorm(20), ncol=2)
X <- cbind(X, 1.5*X[, 1] - pi*X[, 2])
GramSchmidt(X)
Y <- GramSchmidt(X, verbose=TRUE)
round(X - fitted(lm(X ~ Y)), 8)

I'm also thinking about adding a function for the Moore-Penrose inverse as follows (continuing the preceding example):

MoorePenrose <- function(X){
    svd <- SVD(X)
    svd$V %*% diag(1/svd$d) %*% t(svd$U)
}

Y <- MoorePenrose(X)
round(X %*% Y %*% X - X, 8)
round(Y %*% X %*% Y - Y, 8)
round(X %*% Y - t(X %*% Y), 8)
round(Y %*% X - t(Y %*% X), 8)

Thoughts?
John

I'm going to start using matlib shortly in my course, and it would probably save students some grief if I updated the package to CRAN now, to save them having to re-install from github.

Are there any objections to this? Is there anything to do first?

The current version of echelon() just passes the input object to gaussianElimination() with the intention of always expressing the left side as an identity matrix. Should this be replaced with what is commonly meant by the echelon form though? E.g., http://stattrek.com/matrix-algebra/echelon-form.aspx

This relates to my point in issue #11 relating to echelon(..., reduced = FALSE) to distinguish between reduced/non-reduced form.

Hi, I am testing your excellent package but run into a problem! 2D works just fine in Quartz but 3D shows nothing!

test code:

options(rgl.useNULL=TRUE)
library(rgl)
library(matlib)

A <- matrix(c(1, 2, 1, 3, 8, 7, 2, 7, 9), 3, 3)
b <- c(4, 20, 23)

showEqn(A, b)

plotEqn3d(A,b) # Does not work! (Apple Silicon)

Solve(A, b)

There are a few small enhancements since the last release (0.9.3) and one tiny glitch that will soon turn into a nag from Brian Ripley. I've tested the current version with win-builder R-devel (all OK), and also re-built the pkgdown site.

Before I submit, can you take a look and think if there is anything to be added/changed?

matlib 0.9.4

• added xprod() for vector cross-product
• added data-beta vignette
• fix sprintf warning from LU [Phil Chalmers]
• better plotting of planes in plotEqn3d, alpha=0.9

Moving printMatrix and formatNumbers to util.R now generates a warning,

checking for missing documentation entries ... WARNING
Undocumented code objects:
  'formatNumbers' 'printMatrix'

This is because NAMESPACE contains exportPattern("^[[:alpha:]]+"). If they are generally useful enough to be exported, they should be documented.

There has been a function gsorth() in the package, that was copied from heplots for Gram Schmidt orthogonalization, but that is intended for a data frame, and doesn't do what is transparent as GS for a matrix.

I just added the function GramSchmidt() that @philchalmers wrote a while ago. Perhaps it could use a verbose argument.

If this is OK, I will deprecate and/or remove gsorth.

I added the function J() to produce a unit (or constant) vector, matrix, or array. The main use would be in illustrating summation by multiplication by a unit vector on the left or right.

I'm not sure I have this right for the 1D case, and whether or not dimnames= is useful here.

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.

After creating the pkgdown pages, I see that a number of examples now generate warnings. E.g., example(Eigen) gives

....
Eigen> EC$vectors %*% diag(EC$values) %*% t(EC$vectors) # check
    [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    5    6
[3,]    3    6   10
There were 24 warnings (use warnings() to see them)

where the warning is

> warnings()
Warning messages:
1: In (X[, j] %*% E[, k]) * E[, k] :
  Recycling array of length 1 in array-vector arithmetic is deprecated.
  Use c() or as.vector() instead.

These should be tracked down and if possible fixed.

This issue was extracted from #11 so that that long thread can be closed.

The QR function does not seem to be working correctly. It does not return the same value as the base R qr function. Also, the t(q) %*% q doesn't equal I, and q * r doesn't equal A.

Line No. 296
rgl::plot3d(dat, type = 'n', xlim = xlim, ylim = ylim, zlim = c(-3, 3),
xlab = vars[1], ylab = vars[2], zlab = vars[3],
axes=axes)

could/should have 'zlim = zlim' , instead of zlim=c(-3, 3) so user can set its own range for z axis.
default value c(-3, 3) should be set earlier in the code.

I am trying to use the package matlib, however see the type of error I am getting.

If necessary, you can improve the package.

library(matlib) # use the package

A <- matrix(c(0.1, 0.3, 0.5, 0.9, 0.7,0.5), 3, 2)

b <- c(-0.3,-0.1, 0.1)

showEqn(A, b)

plotEqn(A,b)

Solve(A, b) # unique solution

plotEqn(A,b)

0.x[1] + 0.9*x[2] = -0.3

0.3x[1] + 0.7x[2] = -0.1

0.5x[1] + 0.5x[2] = 0.1

Error in parse(text = sub("=", "==", labels[i])) :

:1:5: unexpected symbol

1: 0.x

    ^

Is there a way to re-position the label of the plot lines from the plotEqn() function?

Hi,

I noticed that ?GramSchmidt has this example:

# non-invertible matrix; hence, its basis is not orthonormal
(xx <- matrix(c( 1:3, 3:1, 1, 0, -1), 3, 3))
crossprod(xx)
(zz <- GramSchmidt(xx))
zapsmall(crossprod(zz))

The statement that the basis isn't orthonormal is wrong and contradicted by the last line of the example. Maybe what was meant is that its dimension is less than the number of columns in xx?

John

Dear Michael,

This is nice! Thanks for sharing.

I noticed that you were using webgl in a vignette (linear-equations) and this is something that I want to do so I took a look. It turns out that the webgl is not working in on CRAN https://cran.r-project.org/web/packages/matlib/vignettes/linear-equations.html. I think that CanvasMatrix.js needs to be placed inside inst/doc for this to work. Do you know if CRAN will expose this if the file is put there manually?

Best,

Greg Jefferis.

Currently GramSchmidt contains the line below which removes the zero columns

B <- B[, !apply(B, 2, function(b) all(b == 0))]

however, that makes it hard to relate the result to the input matrix. Suggest that there be an argument to skip over that line so that the output columns directly correspond to the input columns in the same position.

In windows it worked well, but in MacOS, when I load matlib package, there is always an error:

Error : .onLoad failed in loadNamespace() for 'rgl', details:
call: dyn.load(file, DLLpath = DLLpath, ...)
error: unable to load shared object '/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgl/libs/rgl.so':
dlopen(/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgl/libs/rgl.so, 6): Library not loaded: /opt/X11/lib/libGLU.1.dylib
Referenced from: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgl/libs/rgl.so
Reason: image not found
Error: package or namespace load failed for ‘matlib’

In the last CRAN version, we had in DESCRIPTION

Additional_repositories: https://dmurdoch.github.io/drat

which I think was for the rgl examples in vignettes, but I don't recall the details.

devtools::check_rdevel() to winbuilder now flags this with a NOTE

* checking CRAN incoming feasibility ... [12s] NOTE
Maintainer: 'Michael Friendly <[email protected]>'

Availability using Additional_repositories specification:
  ?   ?   https://dmurdoch.github.io/drat

Can this Additional_repositories be removed?

In conjunction with a problem with ugly rendering of 2D ellipses in 3D heplots, Duncan Murdock suggested a solution using a new (to me) rgl property: mask_depth

https://stackoverflow.com/questions/64486375/rgl-how-to-avoid-moire-effect-in-transparent-3d-ellipsoids/64490059

In plotEqn3d, we use alpha=1 by default, so this problem doesn't arise, yet the result of showing regression planes is not very insightful. One can only see the non-occluded portions of the planes.

I've implemented Duncan's suggestion in plotEqn3d() on github, and tested it under different values of alpha. My question is: Should I change the default from alpha=1 to something alpha<1 ? Your thoughts?
[This also affects regec3d where the default results are often uglier]

Extending the example:

> # three consistent equations in three unknowns
> A <- matrix(c(13, -4, 2, -4, 11, -2, 2, -2, 8), 3,3)
> b <- c(1,2,4)
> plotEqn3d(A,b)
> plotEqn(A, b, alpha=0.5)
> plotEqn(A, b, alpha=0.8)

The differences among these are best seen in plots you can rotate.
In static views, I get:

alpha=1

alpha=0.5

alpha=0.8

A simple visualization of how the power method works, at least in the 2D and 3D cases, is a plot of the series of eigenvectors in successive iterations.

But this would require the function to save the history in a returned object.

Solve() prints the result in a nice character form, which it also returns invisibly. As well, it really requires, but doesn't test for, the b argument, so this gives a nonsense result:

> Solve(A)
x1    =  0 
  x2  =  0 
   0  =  1 

For consistency, perhaps our "shadow" functions should return the same form of result as their stats:: counterparts, like Det(), Eigen() do.

This would allow, e.g., Solve(A) %*% b and other uses.

I thought the fix for #27 was trivial. But I had upgraded to R version 3.6.2 and the latest R Studio. This also upgraded Roxygen, required a newer version of RTools, and left me with
(a) all .Rd files re-generated, needing to be committed
(b) errors generated on building vignettes

Error: processing vignette 'gramreg.Rmd' failed with diagnostics: pandoc document conversion failed with error 11 --- failed re-building 'gramreg.Rmd'

I'm not sure how to proceed. I committed my change to the example in GramSchmidt.R to master. But for everything else, I bumped the pkg version to 0.9.3, ran devtools::check() and created this new branch.

Going forward, I think we should try to fix the problems in the devel-0.9.3 branch and then merge into master, but I get errors I don't understand from Travis-CI on this branch.

In the script below, I'm trying to illustrate the idea of projections of a matrix X using two projection matrices, P1 and P2

vals <- c(0, 10)
X <- expand.grid(x1 = vals, x2=vals, x3=vals) |> as.matrix()

# project on just x1, x2 plane
P1 <- rbind(diag(2), c(0,0))
Y1 <- X %*% P1

# oblique projection
P2 <- matrix(c(0.71, 0.71, 0, -0.42, .42, 0.84), ncol=2)
Y2 <- X %*% P2

Part of the illustration is to show the vectors x1, x2, and x3 in (y1, y2) space, e.g.:

Code:

# draw X vectors in Y space
xlim <- ylim <- c(-1.1, 1.1)
axes.x <- c(-1, 1, NA, 0, 0)
axes.y <- c(0, 0, NA, -1, 1)
labs <- c(expression(x[1]), expression(x[2]), expression(x[3]))

op <- par(mar=c(4, 5, 1, 1)+.1)
plot(xlim, ylim, type = "n", asp=1,
     xlab = expression(y[1]), ylab = expression(y[2]),
     cex.lab = 1.8)
plotrix::draw.circle(0, 0, 1, 
                     col = adjustcolor("skyblue", alpha = 0.2))
lines(axes.x, axes.y, col = "grey")
matlib::vectors(P1, labels = labs, cex.lab = 1.8, lwd = 3)

When I run this I get warnings I don't understand or know how to avoid:

1: In arrows(..., angle = ang) :
  zero-length arrow is of indeterminate angle and so skipped
2: In arrows(..., angle = ang) :
  zero-length arrow is of indeterminate angle and so skipped
3: In arrows(..., angle = ang) :
  zero-length arrow is of indeterminate angle and so skipped
4: In arrows(..., angle = ang) :
  zero-length arrow is of indeterminate angle and so skipped
...

As well, I'd like to alter the position of the labels. But, with the following call, no labels are produced:

matlib::vectors(P1, labels = labs, cex.lab = 1.8, lwd = 3, pos.lab = c(4, 3, 0))

Is there a good reason to be using

mass <- requireNamespace("MASS", quietly=TRUE)
if (!mass) stop("fractions=TRUE needs MASS package")
fraction <- MASS::fractions
frac <- function(x) as.character(fraction(x))

instead of just importing MASS::fraction? I was thinking of making a more generic rounding/formatting function to clean up some clutter, but the presence of blocks like this seems a bit strange to me.

Hello,
When although I can install the package, when loading it using library function, I get the following error:
Error in dyn.load(file, DLLpath = DLLpath, ...) :
unable to load shared object '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgl/libs/rgl.so':
dlopen(/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgl/libs/rgl.so, 6): Library not loaded: /opt/X11/lib/libGLU.1.dylib
Referenced from: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgl/libs/rgl.so
Reason: image not found
Error: package or namespace load failed for ‘matlib’:
.onLoad failed in loadNamespace() for 'rgl', details:
call: NULL
error: Loading rgl's DLL failed.
On MacOS, rgl depends on XQuartz, which you can download from xquartz.org.

Coefficient is incorrect in some cases when simplify=TRUE

Reproducible case:
M = matrix(c(-10,2,3,4), ncol=2)
b = c(1,2)
showEqn(M,b)
showEqn(M,b, simplify=TRUE)