Package 'matlib'

Description

These functions are designed mainly for tutorial purposes in teaching & learning matrix algebra ideas and applications to statistical methods using R.

Details

In some cases, functions are provided for concepts available elsewhere in R, but where the function call or name is not obvious. In other cases, functions are provided to show or demonstrate an algorithm, sometimes providing a verbose argument to print the details of computations.

In addition, a collection of functions are provided for drawing vector diagrams in 2D and 3D.

These are not meant for production uses. Other methods are more efficient for larger problems.

Topics

The functions in this package are grouped under the following topics

• Convenience functions:
  tr, R, J, len, vec, Proj, mpower, vandermode

• Determinants: functions for calculating determinants by cofactor expansion
  minor, cofactor, rowMinors, rowCofactors

• Elementary row operations: functions for solving linear equations "manually" by the steps used in row echelon form and Gaussian elimination
  rowadd, rowmult, rowswap

• Linear equations: functions to illustrate linear equations of the form $A x = b$
  showEqn, plotEqn

• Gaussian elimination: functions for illustrating Gaussian elimination for solving systems of linear equations of the form $A x = b$.
  gaussianElimination, Inverse, inv, echelon, Ginv, LU, cholesky, swp

• Eigenvalues: functions to illustrate the algorithms for calculating eigenvalues and eigenvectors
  eigen, SVD, powerMethod, showEig

• Vector diagrams: functions for drawing vector diagrams in 2D and 3D
  arrows3d, corner, arc, pointOnLine, vectors, vectors3d, regvec3d

Most of these ideas and implementations arose in courses and books by the authors. [Psychology 6140](http://friendly.apps01.yorku.ca/psy6140/) was a starting point. Fox (1984) introduced illustrations of vector geometry.

macOS Installation Note

The functions that draw 3D graphs use the rgl package. On macOS, the rgl package requires that XQuartz be installed. After installing XQuartz, it's necessary either to log out of and back into your macOS account or to reboot your Mac.

References

Fox, J. Linear Statistical Models and Related Methods. John Wiley and Sons, 1984

Fox, J. and Friendly, M. (2016). "Visualizing Simultaneous Linear Equations, Geometric Vectors, and Least-Squares Regression with the matlib Package for R". useR Conference, Stanford, CA, June 27 - June 30, 2016.

Description

This function calculates the adjoint of a square matrix, defined as the transposed matrix of cofactors of all elements.

Usage

adjoint(A)

Arguments

Value

a matrix of the same size as A

Author(s)

Michael Friendly

See Also

Other determinants: Det(), cofactor(), minor(), rowCofactors(), rowMinors()

Examples

A <- J(3, 3) + 2*diag(3)
adjoint(A)

Description

angle calculates the angle between two vectors.

Usage

angle(x, y, degree = TRUE)

Arguments

Value

a scalar containing the angle between the vectors

See Also

len

Examples

x <- c(2,1)
y <- c(1,1)
angle(x, y) # degrees
angle(x, y, degree = FALSE) # radians

# visually
xlim <- c(0,2.5)
ylim <- c(0,2)
# proper geometry requires asp=1
plot( xlim, ylim, type="n", xlab="X", ylab="Y", asp=1,
  main = expression(theta == 18.4))
abline(v=0, h=0, col="gray")
vectors(rbind(x,y), col=c("red", "blue"), cex.lab=c(2, 2)) 
text(.5, .37, expression(theta))


####
x <- c(-2,1)
y <- c(1,1)
angle(x, y) # degrees
angle(x, y, degree = FALSE) # radians

# visually
xlim <- c(-2,1.5)
ylim <- c(0,2)
# proper geometry requires asp=1
plot( xlim, ylim, type="n", xlab="X", ylab="Y", asp=1,
  main = expression(theta == 108.4))
abline(v=0, h=0, col="gray")
vectors(rbind(x,y), col=c("red", "blue"), cex.lab=c(2, 2)) 
text(0, .4, expression(theta), cex=1.5)

Description

A utility function for drawing vector diagrams. Draws a circular arc to show the angle between two vectors in 2D or 3D.

Usage

arc(p1, p2, p3, d = 0.1, absolute = TRUE, ...)

Arguments

Details

In this implementation, the two vectors are specified by three points, p1, p2, p3, meaning a line from p1 to p2, and another line from p2 to p3.

Value

none

References

https://math.stackexchange.com/questions/1507248/find-arc-between-two-tips-of-vectors-in-3d

See Also

Other vector diagrams: Proj(), arrows3d(), circle3d(), corner(), plot.regvec3d(), pointOnLine(), regvec3d(), vectors(), vectors3d()

Examples

library(rgl)
vec <- rbind(diag(3), c(1,1,1))
rownames(vec) <- c("X", "Y", "Z", "J")
open3d()
aspect3d("iso")
vectors3d(vec, col=c(rep("black",3), "red"), lwd=2)
# draw the XZ plane, whose equation is Y=0
planes3d(0, 0, 1, 0, col="gray", alpha=0.2)
# show projections of the unit vector J
segments3d(rbind( c(1,1,1), c(1, 1, 0)))
segments3d(rbind( c(0,0,0), c(1, 1, 0)))
segments3d(rbind( c(1,0,0), c(1, 1, 0)))
segments3d(rbind( c(0,1,0), c(1, 1, 0)))
segments3d(rbind( c(1,1,1), c(1, 0, 0)))

# show some orthogonal vectors
p1 <- c(0,0,0)
p2 <- c(1,1,0)
p3 <- c(1,1,1)
p4 <- c(1,0,0)
# show some angles
arc(p1, p2, p3, d=.2)
arc(p4, p1, p2, d=.2)
arc(p3, p1, p2, d=.2)

Description

Draws nice 3D arrows with cone3ds at their tips.

Usage

arrows3d(
  coords,
  headlength = 0.035,
  head = "end",
  scale = NULL,
  radius = NULL,
  ref.length = NULL,
  draw = TRUE,
  ...
)

Arguments

Details

This function is meant to be analogous to arrows, but for 3D plots using rgl. headlength, scale and radius set the length, scale factor and base radius of the arrow head, a 3D cone. The units of these are all in terms of the ranges of the current rgl 3D scene.

Value

invisibly returns the length of the vector used to scale the arrow heads

Author(s)

January Weiner, borrowed from the pca3d package, slightly modified by John Fox

See Also

vectors3d

Other vector diagrams: Proj(), arc(), circle3d(), corner(), plot.regvec3d(), pointOnLine(), regvec3d(), vectors(), vectors3d()

Examples

#none yet

Description

Recover the history of the row operations that have been performed. This function combines the transformation matrices into a single transformation matrix representing all row operations or may optionally print all the individual operations which have been performed.

Usage

buildTmat(x, all = FALSE)

## S3 method for class 'trace'
as.matrix(x, ...)

## S3 method for class 'trace'
print(x, ...)

Arguments

Value

the transformation matrix or a list of individual transformation matrices

Author(s)

Phil Chalmers

See Also

echelon, gaussianElimination

Examples

A <- matrix(c(2, 1, -1,
             -3, -1, 2,
             -2,  1, 2), 3, 3, byrow=TRUE)
b <- c(8, -11, -3)

# using row operations to reduce below diagonal to 0
Abt <- Ab <- cbind(A, b)
Abt <- rowadd(Abt, 1, 2, 3/2)
Abt <- rowadd(Abt, 1, 3, 1)
Abt <- rowadd(Abt, 2, 3, -4)
Abt

# build T matrix and multiply by original form
(T <- buildTmat(Abt))
T %*% Ab    # same as Abt

# print all transformation matrices
buildTmat(Abt, TRUE)

# invert transformation matrix to reverse operations
inv(T) %*% Abt

# gaussian elimination
(soln <- gaussianElimination(A, b))
T <- buildTmat(soln)
inv(T) %*% soln

Description

Returns the Cholesky square root of the non-singular, symmetric matrix X. The purpose is mainly to demonstrate the algorithm used by Kennedy & Gentle (1980).

Usage

cholesky(X, tol = sqrt(.Machine$double.eps))

Arguments

Value

the Cholesky square root of X

Author(s)

John Fox

References

Kennedy W.J. Jr, Gentle J.E. (1980). Statistical Computing. Marcel Dekker.

See Also

chol for the base R function

gsorth for Gram-Schmidt orthogonalization of a data matrix

Examples

C <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric
C
cholesky(C)
cholesky(C) %*% t(cholesky(C))  # check

Description

Draw circles on an existing plot.

Usage

circle(
  x,
  y,
  radius,
  nv = 60,
  border = NULL,
  col = NA,
  lty = 1,
  density = NULL,
  angle = 45,
  lwd = 1
)

Arguments

Details

Rather than depending on the aspect ratio par("asp") set globally or in the call to plot, circle uses the dimensions of the current plot and the x and y coordinates to draw a circle rather than an ellipse. Of course, if you resize the plot the aspect ratio can change.

This function was copied from draw.circle

Value

Invisibly returns a list with the x and y coordinates of the points on the circumference of the last circle displayed.

Author(s)

Jim Lemon, thanks to David Winsemius for the density and angle args

See Also

polygon

Examples

plot(1:5,seq(1,10,length=5),
     type="n",xlab="",ylab="",
     main="Test draw.circle")
# draw three concentric circles
circle(2, 4, c(1, 0.66, 0.33),border="purple",
            col=c("#ff00ff","#ff77ff","#ffccff"),lty=1,lwd=1)
# draw some others
circle(2.5, 8, 0.6,border="red",lty=3,lwd=3)
circle(4, 3, 0.7,border="green",col="yellow",lty=1,
            density=5,angle=30,lwd=10)
circle(3.5, 8, 0.8,border="blue",lty=2,lwd=2)

Description

A utility function for drawing a horizontal circle in the (x,y) plane in a 3D graph

Usage

circle3d(center, radius, segments = 100, fill = FALSE, ...)

Arguments

See Also

Other vector diagrams: Proj(), arc(), arrows3d(), corner(), plot.regvec3d(), pointOnLine(), regvec3d(), vectors(), vectors3d()

Examples

ctr=c(0,0,0)
circle3d(ctr, 3, fill = TRUE)
circle3d(ctr - c(-1,-1,0), 3, col="blue")
circle3d(ctr + c(1,1,0),   3, col="red")

Description

A small artificial data set used to illustrate statistical concepts.

Usage

data("class")

Format

A data frame with 15 observations on the following 4 variables.

sex

a factor with levels F M

age

a numeric vector

height

a numeric vector

weight

a numeric vector

Examples

data(class)
plot(class)

Description

Returns the cofactor of element (i,j) of the square matrix A, i.e., the signed minor of the sub-matrix that results when row i and column j are deleted.

Usage

cofactor(A, i, j)

Arguments

Value

the cofactor of A[i,j]

Author(s)

Michael Friendly

See Also

rowCofactors for all cofactors of a given row

Other determinants: Det(), adjoint(), minor(), rowCofactors(), rowMinors()

Examples

M <- matrix(c(4, -12, -4,
              2,   1,  3,
             -1,  -3,  2), 3, 3, byrow=TRUE)
cofactor(M, 1, 1)
cofactor(M, 1, 2)
cofactor(M, 1, 3)

Description

Draws a cone in 3D from a base point to a tip point, with a given radius at the base. This is used to draw nice arrow heads in arrows3d.

Usage

cone3d(base, tip, radius = 10, col = "grey", scale = NULL, ...)

Arguments

Value

returns the integer object ID of the shape that was added to the scene

Author(s)

January Weiner, borrowed from from the pca3d package

See Also

arrows3d

Examples

# none yet

Description

A utility function for drawing vector diagrams. Draws two line segments to indicate the angle between two vectors, typically used for indicating orthogonal vectors are at right angles in 2D and 3D diagrams.

Usage

corner(p1, p2, p3, d = 0.1, absolute = TRUE, ...)

Arguments

Details

In this implementation, the two vectors are specified by three points, p1, p2, p3, meaning a line from p1 to p2, and another line from p2 to p3.

Value

none

See Also

Other vector diagrams: Proj(), arc(), arrows3d(), circle3d(), plot.regvec3d(), pointOnLine(), regvec3d(), vectors(), vectors3d()

Examples

# none yet

Description

Returns the determinant of a square matrix X, computed either by Gaussian elimination, expansion by cofactors, or as the product of the eigenvalues of the matrix. If the latter, X must be symmetric.

Usage

Det(
  X,
  method = c("elimination", "eigenvalues", "cofactors"),
  verbose = FALSE,
  fractions = FALSE,
  ...
)

Arguments

Value

the determinant of X

Author(s)

John Fox

See Also

det for the base R function

gaussianElimination, Eigen

Other determinants: adjoint(), cofactor(), minor(), rowCofactors(), rowMinors()

Examples

A <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric
A
Det(A)
Det(A, verbose=TRUE, fractions=TRUE)
B <- matrix(1:9, 3, 3) # a singular matrix
B
Det(B)
C <- matrix(c(1, .5, .5, 1), 2, 2) # square, symmetric, nonsingular
Det(C)
Det(C, method="eigenvalues")
Det(C, method="cofactors")

Description

Returns the (reduced) row-echelon form of the matrix A, using gaussianElimination.

Usage

echelon(A, B, reduced = TRUE, ...)

Arguments

Details

When the matrix A is square and non-singular, the reduced row-echelon result will be the identity matrix, while the row-echelon from will be an upper triangle matrix. Otherwise, the result will have some all-zero rows, and the rank of the matrix is the number of not all-zero rows.

Value

the reduced echelon form of X.

Author(s)

John Fox

Examples

A <- matrix(c(2, 1, -1,
             -3, -1, 2,
             -2,  1, 2), 3, 3, byrow=TRUE)
b <- c(8, -11, -3)
echelon(A, b, verbose=TRUE, fractions=TRUE) # reduced row-echelon form
echelon(A, b, reduced=FALSE, verbose=TRUE, fractions=TRUE) # row-echelon form

A <- matrix(c(1,2,3,4,5,6,7,8,10), 3, 3) # a nonsingular matrix
A
echelon(A, reduced=FALSE) # the row-echelon form of A
echelon(A) # the reduced row-echelon form of A

b <- 1:3
echelon(A, b)  # solving the matrix equation Ax = b
echelon(A, diag(3)) # inverting A

B <- matrix(1:9, 3, 3) # a singular matrix
B
echelon(B)
echelon(B, reduced=FALSE)
echelon(B, b)
echelon(B, diag(3))

Description

Eigen calculates the eigenvalues and eigenvectors of a square, symmetric matrix using the iterated QR decomposition

Usage

Eigen(X, tol = sqrt(.Machine$double.eps), max.iter = 100, retain.zeroes = TRUE)

Arguments

Value

a list of two elements: values– eigenvalues, vectors– eigenvectors

Author(s)

John Fox and Georges Monette

See Also

eigen

SVD

Examples

C <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric
C
EC <- Eigen(C) # eigenanalysis of C
EC$vectors %*% diag(EC$values) %*% t(EC$vectors) # check

Description

gaussianElimination demonstrates the algorithm of row reduction used for solving systems of linear equations of the form $A x = B$. Optional arguments verbose and fractions may be used to see how the algorithm works.

Usage

gaussianElimination(
  A,
  B,
  tol = sqrt(.Machine$double.eps),
  verbose = FALSE,
  latex = FALSE,
  fractions = FALSE
)

## S3 method for class 'enhancedMatrix'
print(x, ...)

Arguments

Value

If B is absent, returns the reduced row-echelon form of A. If B is present, returns the reduced row-echelon form of A, with the same operations applied to B.

Author(s)

John Fox

Examples

A <- matrix(c(2, 1, -1,
               -3, -1, 2,
               -2,  1, 2), 3, 3, byrow=TRUE)
  b <- c(8, -11, -3)
  gaussianElimination(A, b)
  gaussianElimination(A, b, verbose=TRUE, fractions=TRUE)
  gaussianElimination(A, b, verbose=TRUE, fractions=TRUE, latex=TRUE)

  # determine whether matrix is solvable
  gaussianElimination(A, numeric(3))

  # find inverse matrix by elimination: A = I -> A^-1 A = A^-1 I -> I = A^-1
  gaussianElimination(A, diag(3))
  inv(A)

  # works for 1-row systems (issue # 30)
  A2 <- matrix(c(1, 1), nrow=1)
  b2 = 2
  gaussianElimination(A2, b2)
  showEqn(A2, b2)
  # plotEqn works for this case
  plotEqn(A2, b2)

Description

Calculate a multiplication factor for the Y dimension to correct for unequal plot aspect and coordinate ratios on the current graphics device.

Usage

getYmult()

Details

getYmult retrieves the plot aspect ratio and the coordinate ratio for the current graphics device, calculates a multiplicative factor to equalize the X and Y dimensions of a plotted graphic object.

Value

The correction factor for the Y dimension.

Author(s)

Jim Lemon

Description

Ginv returns an arbitrary generalized inverse of the matrix A, using gaussianElimination.

Usage

Ginv(A, tol = sqrt(.Machine$double.eps), verbose = FALSE, fractions = FALSE)

Arguments

Details

A generalized inverse is a matrix $\mathbf{A}^-$ satisfying $\mathbf{A A^- A} = \mathbf{A}$.

The purpose of this function is mainly to show how the generalized inverse can be computed using Gaussian elimination.

Value

the generalized inverse of A, expressed as fractions if fractions=TRUE, or rounded

Author(s)

John Fox

See Also

ginv for a more generally usable function

Examples

A <- matrix(c(1,2,3,4,5,6,7,8,10), 3, 3) # a nonsingular matrix
A
Ginv(A, fractions=TRUE)  # a generalized inverse of A = inverse of A
round(Ginv(A) %*% A, 6)  # check

B <- matrix(1:9, 3, 3) # a singular matrix
B
Ginv(B, fractions=TRUE)  # a generalized inverse of B
B %*% Ginv(B) %*% B   # check

Description

Carries out simple Gram-Schmidt orthogonalization of a matrix. Treating the columns of the matrix X in the given order, each successive column after the first is made orthogonal to all previous columns by subtracting their projections on the current column.

Usage

GramSchmidt(
  X,
  normalize = TRUE,
  verbose = FALSE,
  tol = sqrt(.Machine$double.eps),
  omit_zero_columns = TRUE
)

Arguments

Value

A matrix of the same size as X, with orthogonal columns (but with 0 columns removed by default)

Author(s)

Phil Chalmers, John Fox

Examples

(xx <- matrix(c( 1:3, 3:1, 1, 0, -2), 3, 3))
crossprod(xx)
(zz <- GramSchmidt(xx, normalize=FALSE))
zapsmall(crossprod(zz))

# normalized
(zz <- GramSchmidt(xx))
zapsmall(crossprod(zz))

# print steps
GramSchmidt(xx, verbose=TRUE)

# A non-invertible matrix;  hence, it is of deficient rank
(xx <- matrix(c( 1:3, 3:1, 1, 0, -1), 3, 3))
R(xx)
crossprod(xx)
# GramSchmidt finds an orthonormal basis
(zz <- GramSchmidt(xx))
zapsmall(crossprod(zz))

Description

Calculates a matrix with uncorrelated columns using the Gram-Schmidt process

Usage

gsorth(y, order, recenter = TRUE, rescale = TRUE, adjnames = TRUE)

Arguments

Details

This function, originally from the heplots package has now been deprecated in matlib. Use GramSchmidt instead.

Value

a matrix/data frame with uncorrelated columns

Examples

## Not run: 
 set.seed(1234)
 A <- matrix(c(1:60 + rnorm(60)), 20, 3)
 cor(A)
 G <- gsorth(A)
 zapsmall(cor(G))
 
## End(Not run)

Description

Uses gaussianElimination to find the inverse of a square, non-singular matrix, $X$.

Usage

Inverse(X, tol = sqrt(.Machine$double.eps), verbose = FALSE, ...)

Arguments

Details

The method is purely didactic: The identity matrix, $I$, is appended to $X$, giving $X | I$. Applying Gaussian elimination gives $I | X^{-1}$, and the portion corresponding to $X^{-1}$ is returned.

Value

the inverse of X

Author(s)

John Fox

Examples

A <- matrix(c(2, 1, -1,
               -3, -1, 2,
               -2,  1, 2), 3, 3, byrow=TRUE)
  Inverse(A)
  Inverse(A, verbose=TRUE, fractions=TRUE)

Description

This function creates a vector, matrix or array of constants, typically used for the unit vector or unit matrix in matrix expressions.

Usage

J(..., constant = 1, dimnames = NULL)

Arguments

Details

The "dimnames" attribute is optional: if present it is a list with one component for each dimension, either NULL or a character vector of the length given by the element of the "dim" attribute for that dimension. The list can be named, and the list names will be used as names for the dimensions.

Examples

J(3)
J(2,3)
J(2,3,2)
J(2,3, constant=2, dimnames=list(letters[1:2], LETTERS[1:3]))

X <- matrix(1:6, nrow=2, ncol=3)
dimnames(X) <- list(sex=c("M", "F"), day=c("Mon", "Wed", "Fri"))
J(2) %*% X      # column sums
X %*% J(3)      # row sums

Description

len calculates the Euclidean length (also called Euclidean norm) of a vector or the length of each column of a numeric matrix.

Usage

len(X)

Arguments

Value

a scalar or vector containing the length(s)

See Also

norm for more general matrix norms

Examples

len(1:3)
len(matrix(1:9, 3, 3))

# distance between two vectors
len(1:3 - c(1,1,1))

Description

LU computes the LU decomposition of a matrix, $A$, such that $P A = L U$, where $L$ is a lower triangle matrix, $U$ is an upper triangle, and $P$ is a permutation matrix.

Usage

LU(A, b, tol = sqrt(.Machine$double.eps), verbose = FALSE, ...)

Arguments

Details

The LU decomposition is used to solve the equation $A x = b$ by calculating $L(Ux - d) = 0$, where $Ld = b$. If row exchanges are necessary for $A$ then the permutation matrix $P$ will be required to exchange the rows in $A$; otherwise, $P$ will be an identity matrix and the LU equation will be simplified to $A = L U$.

Value

A list of matrix components of the solution, P, L and U. If b is supplied, the vectors $d$ and x are also returned.

Author(s)

Phil Chalmers

Examples

A <- matrix(c(2, 1, -1,
               -3, -1, 2,
               -2,  1, 2), 3, 3, byrow=TRUE)
  b <- c(8, -11, -3)
  (ret <- LU(A)) # P is an identity; no row swapping
  with(ret, L %*% U) # check that A = L * U
  LU(A, b)
  
  LU(A, b, verbose=TRUE)
  LU(A, b, verbose=TRUE, fractions=TRUE)

  # permutations required in this example
  A <- matrix(c(1,  1, -1,
                2,  2,  4,
                1, -1,  1), 3, 3, byrow=TRUE)
  b <- c(1, 2, 9)
  (ret <- LU(A, b))
  with(ret, P %*% A)
  with(ret, L %*% U)

Description

This function provides a soft-wrapper to xtable::xtableMatharray() with additional support for fractions output and brackets.

Usage

matrix2latex(
  x,
  fractions = FALSE,
  brackets = TRUE,
  show.size = FALSE,
  digits = NULL,
  ...
)

Arguments

Details

The code for brackets matches some of the options from the AMS matrix LaTeX package: \pmatrix{}, \bmatrix{}, \Bmatrix{}, ... .

Author(s)

Phil Chalmers

Examples

A <- matrix(c(2, 1, -1,
             -3, -1, 2,
             -2,  1, 2), 3, 3, byrow=TRUE)
b <- c(8, -11, -3)

matrix2latex(cbind(A,b))
matrix2latex(cbind(A,b), digits = 0)
matrix2latex(cbind(A/2,b), fractions = TRUE)

matrix2latex(A, digits=0, brackets="p", show.size = TRUE)

Description

Returns the minor of element (i,j) of the square matrix A, i.e., the determinant of the sub-matrix that results when row i and column j are deleted.

Usage

minor(A, i, j)

Arguments

Value

the minor of A[i,j]

Author(s)

Michael Friendly

See Also

rowMinors for all minors of a given row

Other determinants: Det(), adjoint(), cofactor(), rowCofactors(), rowMinors()

Examples

M <- matrix(c(4, -12, -4,
              2,   1,  3,
             -1,  -3,  2), 3, 3, byrow=TRUE)
minor(M, 1, 1)
minor(M, 1, 2)
minor(M, 1, 3)

Description

The Moore-Penrose inverse is a generalization of the regular inverse of a square, non-singular, symmetric matrix to other cases (rectangular, singular), yet retain similar properties to a regular inverse.

Usage

MoorePenrose(X, tol = sqrt(.Machine$double.eps))

Arguments

Value

The Moore-Penrose inverse of X

Examples

X <- matrix(rnorm(20), ncol=2)
# introduce a linear dependency in X[,3]
X <- cbind(X, 1.5*X[, 1] - pi*X[, 2])

Y <- MoorePenrose(X)
# demonstrate some properties of the M-P inverse
# X Y X = X
round(X %*% Y %*% X - X, 8)
# Y X Y = Y
round(Y %*% X %*% Y - Y, 8)
# X Y = t(X Y)
round(X %*% Y - t(X %*% Y), 8)
# Y X = t(Y X)
round(Y %*% X - t(Y %*% X), 8)

Description

A simple function to demonstrate calculating the power of a square symmetric matrix in terms of its eigenvalues and eigenvectors.

Usage

mpower(A, p, tol = sqrt(.Machine$double.eps))

Arguments

Details

The matrix power p can be a fraction or other non-integer. For example, p=1/2 and p=1/3 give a square-root and cube-root of the matrix.

Negative powers are also allowed. For example, p=-1 gives the inverse and p=-1/2 gives the inverse square-root.

Value

A raised to the power p: A^p

See Also

The {%^%} operator in the expm package is far more efficient

Examples

C <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric
C
mpower(C, 2)
zapsmall(mpower(C, -1))
solve(C)    # check

Description

The plot method for regvec3d objects uses the low-level graphics tools in this package to draw 3D and 3D vector diagrams reflecting the partial and marginal relations of y to x1 and x2 in a bivariate multiple linear regression model, lm(y ~ x1 + x2).

The summary method prints the vectors and their vector lengths, followed by the summary for the model.

Usage

## S3 method for class 'regvec3d'
plot(
  x,
  y,
  dimension = 3,
  col = c("black", "red", "blue", "brown", "lightgray"),
  col.plane = "gray",
  cex.lab = 1.2,
  show.base = 2,
  show.marginal = FALSE,
  show.hplane = TRUE,
  show.angles = TRUE,
  error.sphere = c("none", "e", "y.hat"),
  scale.error.sphere = x$scale,
  level.error.sphere = 0.95,
  grid = FALSE,
  add = FALSE,
  ...
)

## S3 method for class 'regvec3d'
summary(object, ...)

## S3 method for class 'regvec3d'
print(x, ...)

Arguments

Details

A 3D diagram shows the vector y and the plane formed by the predictors, x1 and x2, where all variables are represented in deviation form, so that the intercept need not be included.

A 2D diagram, using the first two columns of the result, can be used to show the projection of the space in the x1, x2 plane.

The drawing functions vectors and link{vectors3d} used by the plot.regvec3d method only work reasonably well if the variables are shown on commensurate scales, i.e., with either scale=TRUE or normalize=TRUE.

Value

None

References

Fox, J. (2016). Applied Regression Analysis and Generalized Linear Models, 3rd ed., Sage, Chapter 10.

See Also

regvec3d, vectors3d, vectors

Other vector diagrams: Proj(), arc(), arrows3d(), circle3d(), corner(), pointOnLine(), regvec3d(), vectors(), vectors3d()

Examples

if (require(carData)) {
   data("Duncan", package="carData")
   dunc.reg <- regvec3d(prestige ~ income + education, data=Duncan)
   plot(dunc.reg)
   plot(dunc.reg, dimension=2)
   plot(dunc.reg, error.sphere="e")
   summary(dunc.reg)

   # Example showing Simpson's paradox
   data("States", package="carData")
   states.vec <- regvec3d(SATM ~ pay + percent, data=States, scale=TRUE)
   plot(states.vec, show.marginal=TRUE)
   plot(states.vec, show.marginal=TRUE, dimension=2)
   summary(states.vec)
}

Description

Shows what matrices $A, b$ look like as the system of linear equations, $A x = b$ with two unknowns, x1, x2, by plotting a line for each equation.

Usage

plotEqn(
  A,
  b,
  vars,
  xlim,
  ylim,
  col = 1:nrow(A),
  lwd = 2,
  lty = 1,
  axes = TRUE,
  labels = TRUE,
  solution = TRUE
)

Arguments

Value

nothing; used for the side effect of making a plot

Author(s)

Michael Friendly

References

Fox, J. and Friendly, M. (2016). "Visualizing Simultaneous Linear Equations, Geometric Vectors, and Least-Squares Regression with the matlib Package for R". useR Conference, Stanford, CA, June 27 - June 30, 2016.

See Also

showEqn

Examples

# consistent equations
A<- matrix(c(1,2,3, -1, 2, 1),3,2)
b <- c(2,1,3)
showEqn(A, b)
plotEqn(A,b)

# inconsistent equations
b <- c(2,1,6)
showEqn(A, b)
plotEqn(A,b)

Description

Shows what matrices $A, b$ look like as the system of linear equations, $A x = b$ with three unknowns, x1, x2, and x3, by plotting a plane for each equation.

Usage

plotEqn3d(
  A,
  b,
  vars,
  xlim = c(-2, 2),
  ylim = c(-2, 2),
  zlim,
  col = 2:(nrow(A) + 1),
  alpha = 0.9,
  labels = FALSE,
  solution = TRUE,
  axes = TRUE,
  lit = FALSE
)

Arguments

Value

nothing; used for the side effect of making a plot

Author(s)

Michael Friendly, John Fox

References

Fox, J. and Friendly, M. (2016). "Visualizing Simultaneous Linear Equations, Geometric Vectors, and Least-Squares Regression with the matlib Package for R". useR Conference, Stanford, CA, June 27 - June 30, 2016.

Examples

# 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)

Description

A utility function for drawing vector diagrams. Find position of an interpolated point along a line from x1 to x2.

Usage

pointOnLine(x1, x2, d, absolute = TRUE)

Arguments

Details

The function takes a step of length d along the line defined by the difference between the two points, x2 - x1. When absolute=FALSE, this step is proportional to the difference, while when absolute=TRUE, the difference is first scaled to unit length so that the step is always of length d. Note that the physical length of a line in different directions in a graph depends on the aspect ratio of the plot axes, and lines of the same length will only appear equal if the aspect ratio is one (asp=1 in 2D, or aspect3d("iso") in 3D).

Value

The interpolated point, a vector of the same length as x1

See Also

Other vector diagrams: Proj(), arc(), arrows3d(), circle3d(), corner(), plot.regvec3d(), regvec3d(), vectors(), vectors3d()

Examples

x1 <- c(0, 0)
x2 <- c(1, 4)
pointOnLine(x1, x2, 0.5)
pointOnLine(x1, x2, 0.5, absolute=FALSE)
pointOnLine(x1, x2, 1.1)

y1 <- c(1, 2, 3)
y2 <- c(3, 2, 1)
pointOnLine(y1, y2, 0.5)
pointOnLine(y1, y2, 0.5, absolute=FALSE)

Description

Finds a dominant eigenvalue, $\lambda_1$, and its corresponding eigenvector, $v_1$, of a square matrix by applying Hotelling's (1933) Power Method with scaling.

Usage

powerMethod(A, v = NULL, eps = 1e-06, maxiter = 100, plot = FALSE)

Arguments

Details

The method is based upon the fact that repeated multiplication of a matrix $A$ by a trial vector $v_1^{(k)}$ converges to the value of the eigenvector,

$v_1^{(k+1)} = A v_1^{(k)} / \vert\vert A v_1^{(k)} \vert\vert$

The corresponding eigenvalue is then found as

$\lambda_1 = \frac{v_1^T A v_1}{v_1^T v_1}$

In pre-computer days, this method could be extended to find subsequent eigenvalue - eigenvector pairs by "deflation", i.e., by applying the method again to the new matrix. $A - \lambda_1 v_1 v_1^{T}$.

This method is still used in some computer-intensive applications with huge matrices where only the dominant eigenvector is required, e.g., the Google Page Rank algorithm.

Value

a list containing the eigenvector (vector), eigenvalue (value), iterations (iter), and iteration history (vector_iterations)

Author(s)

Gaston Sanchez (from matrixkit)

References

Hotelling, H. (1933). Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 24, 417-441, and 498-520.

Examples

A <- cbind(c(7, 3), c(3, 6))
powerMethod(A)
eigen(A)$values[1] # check
eigen(A)$vectors[,1]

# demonstrate how the power method converges to a solution
powerMethod(A, v = c(-.5, 1), plot = TRUE)

B <- cbind(c(1, 2, 0), c(2, 1, 3), c(0, 3, 1))
(rv <- powerMethod(B))

# deflate to find 2nd latent vector
l <- rv$value
v <- c(rv$vector)
B1 <- B - l * outer(v, v)
powerMethod(B1)
eigen(B)$vectors     # check

# a positive, semi-definite matrix, with eigenvalues 12, 6, 0
C <- matrix(c(7, 4, 1,  4, 4, 4,  1, 4, 7), 3, 3)
eigen(C)$vectors
powerMethod(C)

Description

This function is designed to print a collection of matrices, vectors, character strings and matrix expressions side by side. A typical use is to illustrate matrix equations in a compact and comprehensible way.

Usage

printMatEqn(..., space = 1, tol = sqrt(.Machine$double.eps), fractions = FALSE)

Arguments

Value

NULL; A formatted sequence of matrices and matrix operations is printed to the console

Author(s)

Phil Chalmers

See Also

showEqn

Examples

A <- matrix(c(2, 1, -1,
             -3, -1, 2,
             -2,  1, 2), 3, 3, byrow=TRUE)
x <- c(2, 3, -1)

# provide implicit or explicit labels
printMatEqn(AA = A, "*", xx = x, '=', b = A %*% x)
printMatEqn(A, "*", x, '=', b = A %*% x)
printMatEqn(A, "*", x, '=', A %*% x)

# compare with showEqn
b <- c(4, 2, 1)
printMatEqn(A, x=paste0("x", 1:3),"=", b)
showEqn(A, b)

# decimal example
A <- matrix(c(0.5, 1, 3, 0.75, 2.8, 4), nrow = 2)
x <- c(0.5, 3.7, 2.3)
y <- c(0.7, -1.2)
b <- A %*% x - y

printMatEqn(A, "*", x, "-", y, "=", b)
printMatEqn(A, "*", x, "-", y, "=", b, fractions=TRUE)

Description

Print a matrix, allowing fractions or LaTeX output

Usage

printMatrix(
  A,
  parent = TRUE,
  fractions = FALSE,
  latex = FALSE,
  tol = sqrt(.Machine$double.eps)
)

Arguments

Value

The formatted matrix

See Also

fractions

Examples

A <- matrix(1:12, 3, 4) / 6
printMatrix(A, fractions=TRUE)
printMatrix(A, latex=TRUE)

Description

Fitting a linear model, lm(y ~ X), by least squares can be thought of geometrically as the orthogonal projection of y on the column space of X. This function is designed to allow exploration of projections and orthogonality.

Usage

Proj(y, X, list = FALSE)

Arguments

Details

The projection is defined as $P y$ where $P = X (X'X)^- X'$ and $X^-$ is a generalized inverse.

Value

the projection of y on X (if list=FALSE) or a list with elements y and P

Author(s)

Michael Friendly

See Also

Other vector diagrams: arc(), arrows3d(), circle3d(), corner(), plot.regvec3d(), pointOnLine(), regvec3d(), vectors(), vectors3d()

Examples

X <- matrix( c(1, 1, 1, 1, 1, -1, 1, -1), 4,2, byrow=TRUE)
y <- 1:4
Proj(y, X[,1])  # project y on unit vector
Proj(y, X[,2])
Proj(y, X)

# project unit vector on line between two points
y <- c(1,1)
p1 <- c(0,0)
p2 <- c(1,0)
Proj(y, cbind(p1, p2))

# orthogonal complements
y <- 1:4
yp <-Proj(y, X, list=TRUE)
yp$y
P <- yp$P
IP <- diag(4) - P
yc <- c(IP %*% y)
crossprod(yp$y, yc)

# P is idempotent:  P P = P
P %*% P
all.equal(P, P %*% P)

Description

QR computes the QR decomposition of a matrix, $X$, that is an orthonormal matrix, $Q$ and an upper triangular matrix, $R$, such that $X = Q R$.

Usage

QR(X, tol = sqrt(.Machine$double.eps))

Arguments

Details

The QR decomposition plays an important role in many statistical techniques. In particular it can be used to solve the equation $Ax = b$ for given matrix $A$ and vector $b$. The function is included here simply to show the algorithm of Gram-Schmidt orthogonalization. The standard qr function is faster and more accurate.

Value

a list of three elements, consisting of an orthonormal matrix Q, an upper triangular matrix R, and the rank of the matrix X

Author(s)

John Fox and Georges Monette

See Also

qr

Examples

A <- matrix(c(1,2,3,4,5,6,7,8,10), 3, 3) # a square nonsingular matrix
res <- QR(A)
res
q <- res$Q
zapsmall( t(q) %*% q)   # check that q' q = I
r <- res$R
q %*% r                 # check that q r = A

# relation to determinant: det(A) = prod(diag(R))
det(A)
prod(diag(r))

B <- matrix(1:9, 3, 3) # a singular matrix
QR(B)

Description

Returns the rank of a matrix X, using the QR decomposition, QR. Included here as a simple function, because rank does something different and it is not obvious what to use for matrix rank.

Usage

R(X)

Arguments

Value

rank of X

See Also

qr

Examples

M <- outer(1:3, 3:1)
M
R(M)

M <- matrix(1:9, 3, 3)
M
R(M)
# why rank=2?
echelon(M)

set.seed(1234)
M <- matrix(sample(1:9), 3, 3)
M
R(M)

Description

regvec3d calculates the 3D vectors that represent the projection of a two-variable multiple regression model from n-D observation space into the 3D mean-deviation variable space that they span, thus showing the regression of y on x1 and x2 in the model lm(y ~ x1 + x2). The result can be used to draw 2D and 3D vector diagrams accurately reflecting the partial and marginal relations of y to x1 and x2 as vectors in this representation.

Usage

regvec3d(x1, ...)

## S3 method for class 'formula'
regvec3d(
  formula,
  data = NULL,
  which = 1:2,
  name.x1,
  name.x2,
  name.y,
  name.e,
  name.y.hat,
  name.b1.x1,
  name.b2.x2,
  abbreviate = 0,
  ...
)

## Default S3 method:
regvec3d(
  x1,
  x2,
  y,
  scale = FALSE,
  normalize = TRUE,
  name.x1 = deparse(substitute(x1)),
  name.x2 = deparse(substitute(x2)),
  name.y = deparse(substitute(y)),
  name.e = "residuals",
  name.y.hat = paste0(name.y, "hat"),
  name.b1.x1 = paste0("b1", name.x1),
  name.b2.x2 = paste0("b2", name.x2),
  name.y1.hat = paste0(name.y, "hat 1"),
  name.y2.hat = paste0(name.y, "hat 2"),
  ...
)

Arguments

Details

If additional variables are included in the model, e.g., lm(y ~ x1 + x2 + x3 + ...), then y, x1 and x2 are all taken as residuals from their separate linear fits on x3 + ..., thus showing their partial relations net of (or adjusting for) these additional predictors.

A 3D diagram shows the vector y and the plane formed by the predictors, x1 and x2, where all variables are represented in deviation form, so that the intercept need not be included.

A 2D diagram, using the first two columns of the result, can be used to show the projection of the space in the x1, x2 plane.

In these views, the ANOVA representation of the various sums of squares for the regression predictors appears as the lengths of the various vectors. For example, the error sum of squares is the squared length of the e vector, and the regression sum of squares is the squared length of the yhat vector.

The drawing functions vectors and link{vectors3d} used by the plot.regvec3d method only work reasonably well if the variables are shown on commensurate scales, i.e., with either scale=TRUE or normalize=TRUE.

Value

An object of class “regvec3d”, containing the following components

Methods (by class)

• regvec3d(formula): Formula method for regvec3d

• regvec3d(default): Default method for regvec3d

References

Fox, J. (2016). Applied Regression Analysis and Generalized Linear Models, 3rd ed., Sage, Chapter 10.

Fox, J. and Friendly, M. (2016). "Visualizing Simultaneous Linear Equations, Geometric Vectors, and Least-Squares Regression with the matlib Package for R". useR Conference, Stanford, CA, June 27 - June 30, 2016.

See Also

plot.regvec3d

Other vector diagrams: Proj(), arc(), arrows3d(), circle3d(), corner(), plot.regvec3d(), pointOnLine(), vectors(), vectors3d()

Examples

library(rgl)
therapy.vec <- regvec3d(therapy ~ perstest + IE, data=therapy)
therapy.vec
plot(therapy.vec, col.plane="darkgreen")
plot(therapy.vec, dimension=2)

Description

The elementary row operation rowadd adds multiples of one or more rows to other rows of a matrix. This is usually used as a means to solve systems of linear equations, of the form $A x = b$, and rowadd corresponds to adding equals to equals.

Usage

rowadd(x, from, to, mult)

Arguments

Details

The functions rowmult and rowswap complete the basic operations used in reduction to row echelon form and Gaussian elimination. These functions are used for demonstration purposes.

Value

the matrix x, as modified

See Also

echelon, gaussianElimination

Other elementary row operations: rowmult(), rowswap()

Examples

A <- matrix(c(2, 1, -1,
             -3, -1, 2,
             -2,  1, 2), 3, 3, byrow=TRUE)
b <- c(8, -11, -3)

# using row operations to reduce below diagonal to 0
Ab <- cbind(A, b)
(Ab <- rowadd(Ab, 1, 2, 3/2))  # row 2 <- row 2 + 3/2 row 1
(Ab <- rowadd(Ab, 1, 3, 1))    # row 3 <- row 3 + 1 row 1
(Ab <- rowadd(Ab, 2, 3, -4))   # row 3 <- row 3 - 4 row 2
# multiply to make diagonals = 1
(Ab <- rowmult(Ab, 1:3, c(1/2, 2, -1)))
# The matrix is now in triangular form

# Could continue to reduce above diagonal to zero
echelon(A, b, verbose=TRUE, fractions=TRUE)

Description

Returns the vector of cofactors of row i of the square matrix A. The determinant, Det(A), can then be found as M[i,] %*% rowCofactors(M,i) for any row, i.

Usage

rowCofactors(A, i)

Arguments

Value

a vector of the cofactors of A[i,]

Author(s)

Michael Friendly

See Also

Det for the determinant

Other determinants: Det(), adjoint(), cofactor(), minor(), rowMinors()

Examples

M <- matrix(c(4, -12, -4,
              2,   1,  3,
             -1,  -3,  2), 3, 3, byrow=TRUE)
minor(M, 1, 1)
minor(M, 1, 2)
minor(M, 1, 3)
rowCofactors(M, 1)
Det(M)
# expansion by cofactors of row 1
M[1,] %*% rowCofactors(M,1)

Description

Returns the vector of minors of row i of the square matrix A

Usage

rowMinors(A, i)

Arguments

Value

a vector of the minors of A[i,]

Author(s)

Michael Friendly

See Also

Other determinants: Det(), adjoint(), cofactor(), minor(), rowCofactors()

Examples

M <- matrix(c(4, -12, -4,
              2,   1,  3,
             -1,  -3,  2), 3, 3, byrow=TRUE)
minor(M, 1, 1)
minor(M, 1, 2)
minor(M, 1, 3)
rowMinors(M, 1)

Description

Multiplies one or more rows of a matrix by constants. This corresponds to multiplying or dividing equations by constants.

Usage

rowmult(x, row, mult)

Arguments

Value

the matrix x, modified

See Also

echelon, gaussianElimination

Other elementary row operations: rowadd(), rowswap()

Examples

A <- matrix(c(2, 1, -1,
             -3, -1, 2,
             -2,  1, 2), 3, 3, byrow=TRUE)
b <- c(8, -11, -3)

# using row operations to reduce below diagonal to 0
Ab <- cbind(A, b)
(Ab <- rowadd(Ab, 1, 2, 3/2))  # row 2 <- row 2 + 3/2 row 1
(Ab <- rowadd(Ab, 1, 3, 1))    # row 3 <- row 3 + 1 row 1
(Ab <- rowadd(Ab, 2, 3, -4))
# multiply to make diagonals = 1
(Ab <- rowmult(Ab, 1:3, c(1/2, 2, -1)))
# The matrix is now in triangular form

Description

This elementary row operation corresponds to interchanging two equations.

Usage

rowswap(x, from, to)

Arguments

Value

the matrix x, with rows from and to interchanged

See Also

echelon, gaussianElimination

Other elementary row operations: rowadd(), rowmult()

Description

This function is designed for illustrating the eigenvectors associated with the covariance matrix for a given bivariate data set. It draws a data ellipse of the data and adds vectors showing the eigenvectors of the covariance matrix.

Usage

showEig(
  X,
  col.vec = "blue",
  lwd.vec = 3,
  mult = sqrt(qchisq(levels, 2)),
  asp = 1,
  levels = c(0.5, 0.95),
  plot.points = TRUE,
  add = !plot.points,
  ...
)

Arguments

Author(s)

Michael Friendly

See Also

dataEllipse

Examples

x <- rnorm(200)
y <- .5 * x + .5 * rnorm(200)
X <- cbind(x,y)
showEig(X)

# Duncan data
data(Duncan, package="carData")
showEig(Duncan[, 2:3], levels=0.68)
showEig(Duncan[,2:3], levels=0.68, robust=TRUE, add=TRUE, fill=TRUE)

Description

Shows what matrices $A, b$ look like as the system of linear equations, $A x = b$, but written out as a set of equations.

Usage

showEqn(
  A,
  b,
  vars,
  simplify = FALSE,
  reduce = FALSE,
  fractions = FALSE,
  latex = FALSE
)

Arguments

Value

a one-column character matrix, one row for each equation

Author(s)

Michael Friendly, John Fox, and Phil Chalmers

References

Fox, J. and Friendly, M. (2016). "Visualizing Simultaneous Linear Equations, Geometric Vectors, and Least-Squares Regression with the matlib Package for R". useR Conference, Stanford, CA, June 27 - June 30, 2016.

See Also

plotEqn, plotEqn3d

Examples

A <- matrix(c(2, 1, -1,
               -3, -1, 2,
               -2,  1, 2), 3, 3, byrow=TRUE)
  b <- c(8, -11, -3)
  showEqn(A, b)
  # show numerically
  x <- solve(A, b)
  showEqn(A, b, vars=x)

  showEqn(A, b, simplify=TRUE)
  showEqn(A, b, latex=TRUE)

  # lower triangle of equation with zeros omitted (for back solving)
  A <- matrix(c(2, 1, 2,
               -3, -1, 2,
               -2,  1, 2), 3, 3, byrow=TRUE)
  U <- LU(A)$U
  showEqn(U, simplify=TRUE, fractions=TRUE)
  showEqn(U, b, simplify=TRUE, fractions=TRUE)

  ####################
  # Linear models Design Matricies
  data(mtcars)
  ancova <- lm(mpg ~ wt + vs, mtcars)
  summary(ancova)
  showEqn(ancova)
  showEqn(ancova, simplify=TRUE)
  showEqn(ancova, vars=round(coef(ancova),2))
  showEqn(ancova, vars=round(coef(ancova),2), simplify=TRUE)

  twoway_int <- lm(mpg ~ vs * am, mtcars)
  summary(twoway_int)
  car::Anova(twoway_int)
  showEqn(twoway_int)
  showEqn(twoway_int, reduce=TRUE)
  showEqn(twoway_int, reduce=TRUE, simplify=TRUE)

  # Piece-wise linear regression
  x <- c(1:10, 13:22)
  y <- numeric(20)
  y[1:10] <- 20:11 + rnorm(10, 0, 1.5)
  y[11:20] <- seq(11, 15, len=10) + rnorm(10, 0, 1.5)
  plot(x, y, pch = 16)

  x2 <- as.numeric(x > 10)
  mod <- lm(y ~ x + I((x - 10) * x2))
  summary(mod)
  lines(x, fitted(mod))
  showEqn(mod)
  showEqn(mod, vars=round(coef(mod),2))
  showEqn(mod, simplify=TRUE)

Description

Solve the equation system $Ax = b$, given the coefficient matrix $A$ and right-hand side vector $b$, using link{gaussianElimination}. Display the solutions using showEqn.

Usage

Solve(
  A,
  b = rep(0, nrow(A)),
  vars,
  verbose = FALSE,
  simplify = TRUE,
  fractions = FALSE,
  ...
)

Arguments

Details

This function mimics the base function solve when supplied with two arguments, (A, b), but gives a prettier result, as a set of equations for the solution. The call solve(A) with a single argument overloads this, returning the inverse of the matrix A. For that sense, use the function inv instead.

Value

the function is used primarily for its side effect of printing the solution in a readable form, but it invisibly returns the solution as a character vector

Author(s)

John Fox

See Also

gaussianElimination, showEqn inv, solve