CMLS (version 1.0-1)

mlsun: Multivariate Least Squares with Unimodality (and E/I) Constraints

Description

Finds the \(q\) x \(p\) matrix B that minimizes the multivariate least squares problem

sum(( Y - X %*% t(Z %*% B) )^2)

subject to Z %*% B[,j] is unimodal and t(A) %*% B[,j] >= b for all j = 1:p. Unique basis functions and constraints are allowed for each column of B.

Usage

mlsun(X, Y, Z, A, b, meq,
      mode.range = NULL, maxit = 1000, 
      eps = 1e-10, del = 1e-6,
      XtX = NULL, ZtZ = NULL, 
      simplify = TRUE, catchError = FALSE)

Value

If Z is a list with \(q_j = q\) for all \(j = 1,\ldots,p\), then...

B

is returned as a \(q\) x \(p\) matrix when simplify = TRUE

B

is returned as a list of length \(p\) when simplify = FALSE

If Z is a list with \(q_j \neq q\) for some \(j\), then B is returned as a list of length \(p\).

Otherwise B is returned as a \(q\) x \(p\) matrix.

Arguments

X

Matrix of dimension \(n\) x \(p\).

Y

Matrix of dimension \(n\) x \(m\).

Z

Matrix of dimension \(m\) x \(q\). Can also input a list (see Note). If missing, then Z = diag(m) so that \(q = m\).

A

Constraint matrix of dimension \(q\) x \(r\). Can also input a list (see Note). If missing, no equality/inequality (E/I) constraints are imposed.

b

Consraint vector of dimension \(r\) x 1. Can also input a list (see Note). If missing, then b = rep(0, r).

meq

The first meq columns of A are equality constraints, and the remaining r - meq are inequality constraints. Can also input a vector (see Note). If missing, then meq = 0.

mode.range

Mode search ranges, which should be a 2 x \(p\) matrix of integers such that 1 <= mode.range[1,j] <= mode.range[2,j] <= m for all j = 1:p. Default is mode.range = matrix(c(1, m), 2, p).

maxit

Maximum number of iterations for back-fitting algorithm. Ignored if backfit = FALSE.

eps

Convergence tolerance for back-fitting algorithm. Ignored if backfit = FALSE.

del

Stability tolerance for back-fitting algorithm. Ignored if backfit = FALSE.

XtX

Crossproduct matrix: XtX = crossprod(X).

ZtZ

Crossproduct matrix: ZtZ = crossprod(Z).

simplify

If Z is a list, should B be returned as a matrix (if possible)? See Note.

catchError

If catchError = FASLE, an error induced by solve.QP will be returned. Otherwise tryCatch will be used in attempt to catch the error.

Author

Nathaniel E. Helwig <[e-mail]>

Details

A back-fitting algorithm is used to estimate B, where the columns of B are updated sequentially until convergence (outer loop). For each column of B, (the inner loop of) the algorithm searches for the j-th mode across the search range specified by the j-th column of mode.range. The backfitting algorithm is determined to have converged when

mean((B.new - B.old)^2) < eps * (mean(B.old^2) + del),

where B.old and B.new denote the parameter estimates at outer iterations \(t\) and \(t+1\) of the backfitting algorithm.

References

Goldfarb, D., & Idnani, A. (1983). A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Programming, 27, 1-33. tools:::Rd_expr_doi("10.1007/BF02591962")

Helwig, N. E. (in prep). Constrained multivariate least squares in R.

Ten Berge, J. M. F. (1993). Least Squares Optimization in Multivariate Analysis. Volume 25 of M & T Series. DSWO Press, Leiden University. ISBN: 9789066950832

Turlach, B. A., & Weingessel, A. (2019). quadprog: Functions to solve Quadratic Programming Problems. R package version 1.5-8. https://CRAN.R-project.org/package=quadprog

See Also

cmls calls this function for the unimodality constraints.

Examples

######***######   GENERATE DATA   ######***######

# make X
set.seed(2)
n <- 50
m <- 20
p <- 2
Xmat <- matrix(rnorm(n*p), nrow = n, ncol = p)

# make B (which satisfies all constraints except monotonicity)
x <- seq(0, 1, length.out = m)
Bmat <- rbind(sin(2*pi*x), sin(2*pi*x+pi)) / sqrt(4.75)
struc <- rbind(rep(c(TRUE, FALSE), each = m / 2),
               rep(c(FALSE, TRUE), each = m / 2))
Bmat <- Bmat * struc

# make noisy data
set.seed(1)
Ymat <- Xmat %*% Bmat + rnorm(n*m, sd = 0.25)


######***######   UNIMODALITY   ######***######

# unimodal
Bhat.cmls <- cmls(X = Xmat, Y = Ymat, const = "unimod")
Bhat.mlsun <- t(mlsun(X = Xmat, Y = Ymat))
mean((Bhat.cmls - Bhat.mlsun)^2)

# unimodal and structured
Bhat.cmls <- cmls(X = Xmat, Y = Ymat, const = "unimod", struc = struc)
Amat <- vector("list", p)
meq <- rep(0, p)
for(j in 1:p){
   meq[j] <- sum(!struc[j,])
   if(meq[j] > 0){
      A <- matrix(0, nrow = m, ncol = meq[j])
      A[!struc[j,],] <- diag(meq[j])
      Amat[[j]] <- A
   } else {
      Amat[[j]] <- matrix(0, nrow = m, ncol = 1)
   }
}
Bhat.mlsun <- t(mlsun(X = Xmat, Y = Ymat, A = Amat, meq = meq))
mean((Bhat.cmls - Bhat.mlsun)^2)

# unimodal and non-negative
Bhat.cmls <- cmls(X = Xmat, Y = Ymat, const = "uninon")
Bhat.mlsun <- t(mlsun(X = Xmat, Y = Ymat, A = diag(m)))
mean((Bhat.cmls - Bhat.mlsun)^2)

# unimodal and non-negative and structured
Bhat.cmls <- cmls(X = Xmat, Y = Ymat, const = "uninon", struc = struc)
eye <- diag(m)
meq <- rep(0, p)
for(j in 1:p){
   meq[j] <- sum(!struc[j,])
   Amat[[j]] <- eye[,sort(struc[j,], index.return = TRUE)$ix]
}
Bhat.mlsun <- t(mlsun(X = Xmat, Y = Ymat, A = Amat, meq = meq))
mean((Bhat.cmls - Bhat.mlsun)^2)


# see internals of cmls.R for further examples