Nothing
# Copyright 2011-2022, INRAE/INSA/CNRS
# Author: Serguei SOKOL, sokol@insa-toulouse.fr
# License: GPL v2 http://www.gnu.org/licenses/gpl-2.0.html
#' Join elements into a string
#' @description
#' convert elements of vector v (and all following arguments)
#' in strings and join them using sep as separator.
#'
#' @param sep A string used as a separator
#' @param v A string vector to be joined
#' @param ... other variables to be converted to strings and joined
#' @return A joined string
#' @examples
#' join(" ", c("Hello", "World"))
#' @export
join=function(sep, v, ...) {
return(paste(c(v,...), sep="", collapse=sep))
}
#' Non Linear Least Squares with Inequality Constraints
#'
#' @description
#' Solve non linear least squares problem \code{min_par ||r(par,...)$res||}
#' with optional inequality constraints \code{u\%*\%par >= co}
#' and optional equality constraints \code{e\%*\%par = eco}
#' @param par initial values for parameter vector. It can be in non feasible domain,
#' i.e. for which \code{u\%*\%par >= co} is not always respected.
#' @param r a function calculating residual vector
#' by a call to \code{r(par, cjac, ...)}
#' where par is a current parameter vector,
#' and cjac is set to TRUE if Jacobian is required or FALSE if not.
#' A call to \code{r()} must return a list having a field "res" containing the
#' residual vector and optionally a field "jacobian"
#' when cjac is set to TRUE.
#' @param ... parameters passed through to r()
#' @param u constraint matrix in \code{u\%*\%par >= co}
#' @param co constraint right hand side vector
#' @param control a list of control parameters ('=' indicates default values):
#' - errx=1.e-7 error on l2 norm of the iteration step sqrt(pt*p).
#' - maxit=100 maximum of newton iterations
#' - btstart=1 staring value for backtracking coefficient
#' - btfrac=0.5 (0;1) by this value we diminish the step till tight up
#' to the quadratic model of norm reduction in backtrack (bt) iterations
#' - btdesc=0.1 (0;1) how good we have to tight up to the quadratic model
#' 0 - we are very relaxe, 1 - we are very tight (may need many bt iterations)
#' - btmaxit=15 maximum of backtrack iterations
#' - btkmin=1.e-7 a floor value for backtracking fractioning
#' - trace=0 no information is printed during iterations (1 - print is active)
#' - rcond=1.e10 maximal condition number in QR decomposition
#' starting from which a matrix is considered as numerically rank
#' deficient. The inverse of this number is also used as a measure of very
#' small number.
#' - ci = list of options relative to confidence interval estimation
#' + p=0.95 p-value of confidence interval. If you need a plain sd value,
#' set p-value to 0.6826895
#' + report=T report (or not and hence calculate or not) confidence
#' interval information (in the field hci, as 'half confidence interval' width)
#' - history=TRUE report or not the convergence history
#' - adaptbt=FALSE use or not adaptive backtracking
#' - mnorm=NULL a norm matrix for a case sln==TRUE (cf next parameter)
#' - sln=FALSE use or not (default) least norm of the solution
#' (instead of least norm of the increment)
#' - maxstep=NULL a maximal norm for increment step (if not NULL), must be positive
#' - monotone=FALSE should or not the cost decrease be monotone. If TRUE, then at
#' first non decrease of the cost, the iterations are stopped with a warning message.
#' @param e linear equality constraint matrix in \code{e\%*\%par = eco} (dense)
#' @param eco right hand side vector in equality constraints
#' @param flsi function solving linear least squares problem. For a custom function,
#' see interfaces in \code{lsi} or \code{lsi_ln} help pages.
#' @return a list with following components (some components can be absent depending
#' on 'control' parameter)
#' - 'par' estimated values of par as vector
#' - 'paro' the same but in original structure (i.e. matrix if par is a matrix)
#' - 'lastp' the last LSI solution during non linear iterations
#' - 'hci' vector of half-width confidence intervals for par
#' - 'ci_p' p-value for which CI was calculated
#' - 'ci_fdeg' freedom degree used for CI calculation
#' - 'sd_res' standard deviation of residuals
#' - 'covpar' covariance matrix for par
#' - 'laststep' the last increment after possible back-tracking iterations
#' - 'normp' the norm of lastp
#' - 'res' the last residual vector
#' - 'prevres' residual vector from previous non linear iteration
#' - 'jacobian' the last used jacobian
#' - 'retres' last returned result of \code{r()} call
#' - 'it' non linear iteration number where solution was obtained
#' - 'btit' back-tracking iteration number done during the last non linear iteration
#' - 'history' list with convergence history information
#' - 'error' error code: 0 - normal end, 1 - some error occurred, see message in 'mes'
#' - 'mes' textual message explaining what problem was in case of error
#' @details
#' Solving method consist in sequential LSI problems globalized by backtracking technique.
#' If e, eco are not NULL, reduce jacobian to basis of e's kernel before \code{lsi()} call.\cr
#' NB. If the function \code{r()} returns a list having a field "jacobian" it is
#' supposed to be equal to the jacobian dr/dpar.
#' If not, numerical derivation numDeriv::jacobian()
#' is automatically used for its calculation.\cr
#' NB2. nlsic() does not call stop() on possible errors. Instead, 'error' field is set to 1
#' in the returned result. This is done to allow a user to examine the current state
#' of data and possibly take another path then to merely stop the program. So, a
#' user must allways check this field at return from nlsic().\cr
#' NB3. User should test that field 'mes' is not NULL even when error is 0. It may
#' contain a warning message.
#' @seealso [lsi], [lsi_ln], [uplo2uco]
#' @examples
#' # solve min_{a,b} ||exp(a*x+b)-meas||, a,b>=1
#' a_true=1; b_true=2; x=0:5
#' # simulation function
#' sim=function(par, x) exp(par[["a"]]*x+par[["b"]])
#' # residual function to be passed to nlsic()
#' resid=function(par, cjac, ...) {
#' dots=list(...)
#' s=sim(par, dots$x)
#' result=list(res=s-dots$meas)
#' if (cjac) {
#' result$jacobian=cbind(a=s*dots$x, b=s)
#' }
#' result
#' }
#' # simulated noised measurements for true parameters
#' set.seed(7) # for reproducible results
#' meas=sim(c(a=a_true, b=b_true), x)+rnorm(x)
#' # starting values for par
#' par=c(a=0, b=0)
#' # prepare constraints
#' uco=uplo2uco(par, lower=c(a=1, b=1))
#' # main call: solve the problem
#' fit=nlsic(par, resid, uco$u, uco$co, control=list(trace=1), x=x, meas=meas)
#' if (fit$error == 1) {
#' stop(fit$mes)
#' } else {
#' print(fit$par) # a=1.001590, b=1.991194
#' if (! is.null(fit$mes)) {
#' warning(fit$mes)
#' }
#' }
#' @export
nlsic=function(par, r, u=NULL, co=NULL, control=list(), e=NULL, eco=NULL, flsi=lsi, ...) {
#print(match.call())
n=length(par)
m=NROW(co)
nm_par=names(par)
# predefined controls, then overwritten by actual values
con=list(errx=1.e-7, maxit=100, btstart=1., btfrac=0.5, btdesc=0.1,
btmaxit=15, btkmin=1.e-7, trace=0, rcond=1.e10, ci=list(p=0.95, report=T),
history=FALSE, adaptbt=FALSE, mnorm=NULL, sln=FALSE, maxstep=NULL, monotone=FALSE)
nmsC=names(con)
nmsci=names(con$ci)
# predefined confidence interval parameters, then overwritten
ci=con$ci
if (!is.null(control$ci)) {
ci[(namci <- names(control$ci))] <- control$ci
} else {
namci <- names(con$ci)
}
con[(namc <- names(control))] <- control
con$ci=ci
tol=1.e-10
mes=NULL
if (length(noNms <- namc[!namc %in% nmsC])) {
mes=join("", "nlsic: unknown name(s) in control: ", join(", ", noNms))
return(list(par=NULL, error=1, mes=mes))
}
if (length(noNms <- namci[!namci %in% nmsci])) {
mes=join("", "nlsic: unknown names in control$ci: ", join(", ", noNms))
return(list(par=NULL, error=1, mes=mes))
}
if (!is.null(con$maxstep) && con$maxstep <= 0.) {
mes=sprintf("nlsic: maxstep parameter in control list must be positive, got %g instead", con$maxstep)
return(list(par=NULL, error=1, mes=mes))
}
if (!is.null(con$maxstep) && con$maxstep <= con$errx) {
mes=sprintf("nlsic: maxstep parameter in control list is less than errx (maxstep=%g, errx=%g)", con$maxstep, con$errx)
return(list(par=NULL, error=1, mes=mes))
}
# history matrices
hpar=hp=hstep=hres=c()
if (con$history) {
hpar=cbind(received=par)
hist=list()
} else {
hist=NULL
}
if (!is.null(e)) {
e=matrix(e, ncol=n)
ne=nrow(e)
} else {
ne=0
}
econstr=!is.null(e) && !is.null(eco) && nrow(e)>0 && ncol(e)==length(par)
# make u matrix if u=NULL
if (is.null(u)) {
u=matrix(0, nrow=0, ncol=n)
co=c()
}
if (econstr) {
nr=nrow(e)
nc=ncol(e)
if (nr > nc) {
return(list(par=NULL, error=1, mes="nlsic: matrix e is over determined."))
}
# affine transform epar -> par s.t. e*par=eco
# par=parp+Null(t(e))%*%epar, where epar is new parameter vector
nte=Nulla(t(e))
qe=attr(nte, "qr")
if (qe$rank < nr) {
return(list(par=NULL, error=1, mes="nlsic: matrix e is rank deficient."))
}
# particular solution
parp=qr.qy(qe, c(backsolve(qe$qr, eco[qe$pivot], transpose=T), double(n-nr)))
ue=u%*%nte
epar=crossprod(nte, c(par)-parp)
par[]=nte%*%epar+parp
}
ine=if (nrow(u) > 0) co-u%*%c(par) else NULL
if (!is.null(ine) && any(ine > 1./con$rcond)) {
# project par into feasible domain
# in case if residual function is defined only on feasible domain
if (econstr) {
x=ldp(ue, ine, con$rcond)
if (is.null(x)) {
return(list(
par=NULL,
error=1,
mes="nlsic: Unfeasible constraints at starting point"
))
}
laststep=nte%*%x
par[]=par+c(laststep)
} else {
x=ldp(u, ine, con$rcond)
if (is.null(x)) {
mes="nlsic: Unfeasible constraints at starting point (u%*%par-co>0):"
if (!is.null(rownames(u))) {
mes=join("\n", c(mes, paste(rownames(u), " (", -ine, ")", sep="")))
}
return(list(
par=NULL,
error=1,
mes=mes
))
}
laststep=x
par[]=par+c(laststep)
}
}
# newton globalized iterations
it=0
btit=0
converged=F
laststep=par
laststep[]=0.
p=laststep
a=res=lres=vres=normp=NULL
while (!converged && it < con$maxit) {
mes=NULL
if (it == 0) {
#browser()
# residual vector
#cat("jac it=0", par, sep="\n")
lres=r(par, cjac=TRUE, ...)
if (!is.null(lres$err) && lres$err) {
return(list(
par=c(par),
paro=par,
retres=lres,
it=it,
btit=btit,
error=1,
mes=paste("nlsic: Problem in first residual calculation",
lres$mes, sep="\n"))
)
}
res=as.numeric(lres$res)
iva=which(!is.na(res)) # valid entries in res
if (length(iva)==0) {
# we check it only at it=0 hoping that NA-structure remains the same
# all the long ietrations
return(list(
par=NULL,
error=1,
mes="nlsic: no valid residual value."
))
}
vres=res[iva]
b=-res
norb=sqrt(norm2(vres))
#cat("norb0=", norb, "\n", sep="")
norres=norb
klast=0
if (con$trace) {
cat("it=", it, "\tres=", norb, "\n", sep="")
}
if (norb < 1.e-20) {
# we are already at solution point
converged=TRUE
normp=0
a=NULL
next
}
} else {
b=-res
norb=norres
}
# jacobian
a=NULL
if ("jacobian" %in% names(lres) && !is.null(lres$jacobian)) {
a=lres$jacobian[iva,,drop=FALSE]
colnames(a)=nm_par
} else {
# may be we are here because the last call was with cjac=FALSE
if (it > 0) {
#cat("jac", par, sep="\n")
lres=r(par, cjac=TRUE, ...)
res=as.numeric(lres$res)
iva=!is.na(res)
vres=res[iva]
b=-res
norres=sqrt(norm2(vres))
norb=norres
if ("jacobian" %in% names(lres) && !is.null(lres$jacobian)) {
a=lres$jacobian[iva,,drop=FALSE]
colnames(a)=nm_par
}
}
}
if (is.null(a)) {
if (requireNamespace("numDeriv", quietly = TRUE)) {
a=numDeriv::jacobian(function(x) as.numeric(r(x, cjac=FALSE, ...)$res), par)[iva,,drop=FALSE]
colnames(a)=nm_par
} else {
return(list(
par=c(par),
paro=par,
it=it,
btit=btit,
error=1,
history=hist,
jacobian=a,
retres=lres,
mes=paste("nlsic: Not found 'jacobian' in results of 'r()' when cjac=TRUE and numDeriv package is not available.",
attr(p, "mes"), sep="\n")))
}
}
# solve linear ls
#browser()
ine=if (nrow(u) > 0) co-u%*%c(par) else NULL
# newton direction
if (econstr) {
if (con$sln) {
p=flsi(a%*%nte, -vres, ue, ine, rcond=con$rcond, con$mnorm, -c(par))
} else {
p=flsi(a%*%nte, -vres, ue, ine, rcond=con$rcond)
}
p=nte%*%p
} else {
if (con$sln) {
p=flsi(a, -vres, u, ine, rcond=con$rcond, con$mnorm, -c(par))
} else {
p=flsi(a, -vres, u, ine, rcond=con$rcond)
}
}
if (anyNA(p)) {
names(par)=nm_par
return(list(
par=c(par),
paro=par,
it=it,
btit=btit,
error=1,
history=hist,
jacobian=a,
retres=lres,
step=p,
mes=paste("nlsic: NA found in the step",
attr(p, "mes"), sep="\n")))
} else if (! is.null(attr(p, "mes"))) {
mes=join("\n", mes, attr(p, "mes"))
}
if (!all(u%*%(c(par)+c(p))-co >= -tol)) {
names(par)=nm_par
return(list(
par=c(par),
paro=par,
it=it,
btit=btit,
error=1,
history=hist,
jacobian=a,
retres=lres,
step=p,
mes=paste("nlsic: found step does not satisfy inequalities",
attr(p, "mes"), sep="\n")))
}
names(p)=nm_par
#browser()
if (con$history) {
hpar=cbind(hpar, par)
colnames(hpar)[ncol(hpar)]=it
hp=cbind(hp, c(p))
hres=cbind(hres, res)
hist=list(par=hpar, p=hp, res=hres, step=hstep)
}
normp=sqrt(norm2(p))
converged=(normp <= con$errx)
if (normp == 0) {
############### no need for backtracking at this stage
laststep=p
par[]=par+c(laststep)
it=it+1
btit=0
resdecr=NULL
#cat("normp=0", par, sep="\n")
res=as.numeric(r(par, cjac=FALSE, ...)$res)
iva=!is.na(res)
vres=res[iva]
norres=sqrt(norm2(vres))
if (con$history) {
hstep=cbind(hstep, c(laststep))
hist=list(par=hpar, p=hp, res=hres, step=hstep)
}
if (con$trace && (it < 10 || !it%%10)) {
cat("it=", it, "\tres=", norres, "\tnormstep=", normp, "\tbtk=", klast, "\n", sep="")
}
break
}
# check if the direction p is descending
ap=as.double(a%*%p)
n2ap=norm2(ap)
resap=crossprod(vres, ap)
if (resap > 0.) {
laststep=p*NA
mes=paste("nlsic: LSI returned not descending direction",
attr(p, "mes"), sep="\n")
break
}
k=con$btstart; # fraction of p
if (!is.null(con$maxstep)) {
kms=con$maxstep/normp
if (kms < k) {
k=kms
}
}
# backtrack iterations
btit=0
descending=FALSE
while (!descending && btit < con$btmaxit && k >= con$btkmin) {
laststep=k*p
#cat("btit=", btit, par+c(laststep), sep="\n")
lres=try(r(par+c(laststep), cjac=FALSE, ...))
if (inherits(lres, "try-error")) {
k=max(k*con$btfrac, con$btkmin)
btit=btit+1
next
}
if (isTRUE(lres$err != 0)) {
return(list(
par=c(par),
laststep=c(laststep),
retres=lres,
it=it,
btit=btit,
error=1,
mes=paste("nlsic: Problem in residual calculation",
lres$mes, sep="\n"))
)
}
res=as.numeric(lres$res)
iva=!is.na(res) # valid entries in res
vres=res[iva]
norres=sqrt(norm2(vres))
#cat("norres=", norres, "\n", sep="")
#scaprod=crossprod(vres, ap)-resap
#resdecr=sqrt(scaprod); # residual decreasing
realdecr=crossprod(b[iva]-vres, vres+b[iva]) # normally it is positive
lindecr=-k*(2*resap+k*n2ap)
#descending=(
# scaprod >=0. &&
# sqrt(scaprod) >= con$btdesc*sqrt(n2ap*k) &&
# norres < norb
#)
descending=(realdecr>=con$btdesc*lindecr)
klast=k
if (con$adaptbt) {
h=max(min(2., 1./(2.-(norres-norb)*(norres+norb)/(resap*k))), con$btfrac)
# check that p*k*h won't violate inequalities (for h > 1 only)
if (h > 1. && nrow(u) > 0) {
if (any(u%*%(c(par)+c((k*h)*p))-co <= -tol)) {
h=con$btfrac
}
}
k=k*h
#cat("h=", h, "klast=", klast, "k=", k, "\n")
} else {
k=k*con$btfrac
}
k=max(k, con$btkmin)
btit=btit+1
}
if (inherits(lres, "try-error")) {
return(list(
par=c(par),
paro=par,
laststep=rep(NA, length(par)),
retres=lres,
it=it,
btit=btit,
error=1,
mes="nlsic: residual calculation failed during all backtracking iterations"
))
}
par[]=par+c(laststep)
it=it+1
if (con$monotone) {
if (norres > norb) {
converged=T
mes=join("\n", mes, "nlsic: Non monotone descrease of cost function has occured.")
}
}
if (con$trace && ((it < 10 || !it%%10) || converged)) {
cat("it=", it, "\tres=", norres, "\tnormstep=", normp, "\tbtk=", klast, "\n", sep="")
}
if (con$history) {
hstep=cbind(hstep, laststep)
hist=list(par=hpar, p=hp, res=hres, step=hstep)
}
}
if (con$trace && !(it < 10 || !it%%10) && !converged) {
cat("it=", it, "\tres=", norres, "\tnormstep=", normp, "\tbtk=", klast, "\n", sep="")
}
if (it >= con$maxit) {
mes=join("\n", mes, "nlsic: Maximal non linear iteration number is achieved")
}
if (btit >= con$btmaxit) {
mes=join("\n", mes, "nlsic: Maximal backtrack iteration number is achieved")
}
# restore names
names(par)=names(laststep)=names(p)=nm_par
# confidence interval
#browser()
nr=length(vres)
fdeg=max(1., nr-n+ne) # freedom degrees
if (con$ci$report) {
sd_res=sqrt(sum(vres**2)/fdeg)
if (is.null(a)) {
#cat("first jac", par, sep="\n")
a=r(par, TRUE, ...)$jacobian
if (is.null(a)) {
#cat("num jac", par, sep="\n")
if (requireNamespace("numDeriv", quietly = TRUE)) {
a=numDeriv::jacobian(function(x) as.numeric(r(x, cjac=FALSE, ...)$res), par)
colnames(a)=nm_par
} else {
return(list(
par=c(par),
paro=par,
it=it,
btit=btit,
error=1,
history=hist,
jacobian=a,
retres=lres,
mes=paste("nlsic: Not found 'jacobian' in results of 'r()' when cjac=TRUE and numDeriv package is not available.",
attr(p, "mes"), sep="\n")))
}
}
}
if (econstr) {
jac=a%*%nte
} else {
jac=a
}
covpar=try(solve(crossprod(jac)), silent=T)
if (inherits(covpar, "try-error")) {
sv=svd(jac)
covpar=tcrossprod(sv$v*rep(1./pmax(sv$d, 1./con$rcond), rep(nrow(sv$v), ncol(sv$v))))
}
if (econstr) {
covpar=nte%*%tcrossprod(covpar, nte)
}
dimnames(covpar)=list(nm_par, nm_par)
d=sqrt(diag(covpar))
hci=(-stats::qt((1.-con$ci$p)/2., fdeg)*sd_res)*d
names(hci)=nm_par
} else {
sd_res=hci=covpar=ci_fdeg=NULL
}
return(list(
par=c(par),
paro=par,
lastp=c(p),
hci=hci,
ci_p=con$ci$p,
ci_fdeg=fdeg,
sd_res=sd_res,
covpar=covpar,
laststep=c(laststep),
normp=normp,
res=c(res),
prevres=c(vres),
jacobian=a,
retres=lres,
it=it,
btit=btit,
history=hist,
error=0,
mes=mes)
)
}
#' Linear Least Squares with Inequality constraints (LSI)
#'
#' @description
#' solve linear least square problem (min ||A%*%x-b||)
#' with inequality constraints \code{u%*%x>=co}
#' @param a dense matrix A or its QR decomposition
#' @param b right hand side vector. Rows containing NA are dropped.
#' @param u dense matrix of inequality constraints
#' @param co right hand side vector of inequality constraints
#' @param rcond maximal condition number for determining rank deficient matrix
#' @param mnorm dummy parameter
#' @param x0 dummy parameter
#' @return solution vector whose attribute 'mes' may contain a message about possible numerical problems
#' @seealso [lsi_ln], [ldp], [base::qr]
#' @details
#' Method:\cr
#' 1. reduce the problem to ldp (min(xat*xa) => least distance programming)\cr
#' 2. solve ldp\cr
#' 3. change back to x\cr
#' If b is all NA, then a vector of NA is returned.
#'
#' mnrom, and x0 are dummy parameters which are here to make lsi()
#' compatible with lsi_ln() argument list
#' @export
lsi=function(a, b, u=NULL, co=NULL, rcond=1.e10, mnorm=NULL, x0=NULL) {
tol=1.e-10
# remove NA from b
i0=which(is.na(b))
if (length(i0)) {
if (length(i0) == length(b)) # all NA
return(rep(NA_real_, if (is.qr(a)) ncol(a$qr) else ncol(a)))
if (is.qr(a))
a=qr.X(a)
a=a[-i0,,drop=FALSE]
b=b[-i0]
}
if (! is.qr(a)) {
n=ncol(a)
aq=base::qr(a, LAPACK=TRUE)
} else {
n=ncol(a$qr)
aq=a
}
d=abs(diag(aq$qr))
#cat("d=", d, "\n")
aq$rank=sum(d>d[1]/rcond)
if (is.na(aq$rank)) {
x=rep(NA, n)
attr(x, "mes")="lsi: Rank could not be estimated in least squares\n"
return(x)
}
if (aq$rank < n) {
x=rep(NA, n)
attr(x, "mes")=
paste("lsi: Rank deficient matrix in least squares\n",
n-aq$rank,
" unsolvable variable(s):\n",
paste(dimnames(a)[[2]][aq$pivot[-(1:aq$rank)]],
aq$pivot[-(1:aq$rank)], sep="\t", collapse="\n"), "\n",
sep="")
return(x)
}
#x0=qr.solve(aq, b)
x0=qr.coef(aq, b)
if (!is.null(u) && nrow(u)>0) {
# we do have inequalities
# prepare variable change
ut=t(base::backsolve(aq$qr, t(as.matrix(u[, aq$pivot, drop=FALSE])), transpose=TRUE))
xa=ldp(ut, co-u%*%x0, rcond)
if (is.null(xa)) {
# Infeasible constraints detected by ldp
xa=rep(NA, n)
attr(xa, "mes")="lsi: ldp() revealed unfeasible constrants"
return(xa)
}
x=xa
xa=base::backsolve(aq$qr, xa)
x[aq$pivot]=xa
x=x+x0
cou=co-u%*%x
if (any(cou > tol)) {
# round errors made fail inequality constraints
# force them as much as possible even by degrading the residual
dx=ldp(u, cou, rcond)
if (!is.null(dx)) {
x=x+dx
} # else leave x as is
}
#stopifnot(all(u%*%x-co >= -tol))
} else {
# plain QR
x=x0
}
return(x)
}
#' Linear Least Squares with Inequality constraints, least norm solution
#'
#' @description
#' solve linear least square problem \code{min_x ||A*x-b||}
#' with inequality constraints \code{u\%*\%x >= co}
#' If A is rank deficient, least norm solution \code{||mnorm\%*\%(x-x0)||} is used.
#' If the parameter mnorm is NULL, it is treated as an identity matrix.
#' If the vector x0 is NULL, it is treated as 0 vector.
#' @param a dense matrix A or its QR decomposition
#' @param b right hand side vector
#' @param u dense matrix of inequality constraints
#' @param co right hand side vector of inequality constraints
#' @param rcond maximal condition number for determining rank deficient matrix
#' @param mnorm norm matrix (can be dense or sparse) for which \code{\%*\%} operation with a dense vector is defined
#' @param x0 optional vector from which a least norm distance is searched for
#' @return solution vector whose attribute 'mes' may contain a message about possible numerical problems
#' @seealso [lsi], [ldp], [base::qr]
#' @export
lsi_ln=function(a, b, u=NULL, co=NULL, rcond=1.e10, mnorm=NULL, x0=NULL) {
#if(anyNA(u))
# browser()
mes=""
tol=1.e-10
if (! is.qr(a)) {
a=as.matrix(a)
n=ncol(a)
aq=base::qr(a, LAPACK=TRUE)
} else {
n=ncol(a$qr)
aq=a
}
d=diag(aq$qr)
aq$rank=sum(abs(d)>abs(d[1])/rcond)
if (aq$rank == 0) {
# matrix a is zero => no x can diminish the residual norm
if (!is.null(u) && nrow(u) > 0) {
# there are some inequalities to deal with
if (any(co > tol)) {
# at least we satisfy inequalities
x=ldp(u, co, rcond)
if (is.null(x)) {
x=rep(NA, n)
attr(x, "mes")="lsi_ln: matrix a is zero rank and ldp() revealed unfeasible constrants"
}
#stopifnot(all(u%*%x-co >= -tol))
return(x)
} else {
# they are already satisfied => return 0 vector
return(double(n))
}
} else {
# no inequalities
return(double(n))
}
}
# here after the rank is > 0
rdefic=aq$rank < n
if (!rdefic) {
# just passthrough params to lsi()
return(lsi(aq, b, u, co, rcond))
}
# prepare free variable substitution
ndefic=n-aq$rank
i1=seq(len=aq$rank)
ndefic=n-aq$rank
i2=seq(aq$rank+1, n, len=ndefic)
r=base::qr.R(aq)[i1,,drop=FALSE]
r1=r[,i1,drop=FALSE]
r2=r[,i2,drop=FALSE]
qrr=base::qr(t(r), LAPACK=T)
# null space and its complement bases
ba=base::qr.Q(qrr, complete=T)
bal=ba[,i1,drop=FALSE]
ban=ba[,i2,drop=FALSE]
# least norm without constraints
x=bal%*%base::backsolve(qrr$qr, qr.qty(aq, b)[qrr$pivot], transpose=TRUE)
#x0=backsolve(r1, bt) # implicitly, free variables are set to zero
#x=c(x0, rep(0., ndefic)) # we unpivot just before return
#browser()
if (!is.null(u) && nrow(u)>0) {
up=u[,aq$pivot,drop=FALSE]
cou=co-up%*%x
if (any(cou > tol)) {
uz=up%*%ban
# minimize ||y;z||
ut=t(base::backsolve(qrr$qr, t(as.matrix(up%*%bal)), transpose=FALSE))
ut=cbind(ut, uz)
ya=ldp(ut, cou, rcond)
if (!is.null(ya)) {
# move along the borders to diminish ||a*x-b||
# the solution is searched for in the form
# y=ya+bau*yz
# ia are indexes of border equations to add as equality constraints (active set)
cou=as.numeric(cou-ut%*%ya)
ia=which(cou > -tol) # must not be empty if we are here
if (length(ia) == 0L) {
# but it can happen for badly conditioned ut
# let take the inequality closest to 0
ia=which.max(cou)
}
bau=Nulla(t(ut[ia,,drop=FALSE]))
#if (anyNA(ut[-ia,,drop=FALSE]%*%bau))
# cat("sent NA in lsi_ln(.., u, ...)\n")
yz=lsi_ln(bau[i1,,drop=FALSE], -ya[i1], ut[-ia,,drop=FALSE]%*%bau, cou[-ia], rcond=rcond)
if (!anyNA(yz)) {
# got a solution
ya=ya+bau%*%yz
dx=bal%*%base::backsolve(qrr$qr, ya[i1], transpose=T)+ban%*%ya[i2]
x=x+dx
attr(x, "mes")=attr(yz, "mes")
} else {
x=rep(NA, n)
attr(x, "mes")=join("\n", attr(x, "mes"), "lsi_ln: failed to glide along borders.")
return(x)
}
} else {
# unfeasible constraints
x=rep(NA, n)
attr(x, "mes")="lsi_ln: unfeasible constraints"
return(x)
}
# solve lsi_ln within just Null space (without modifying ||y||)
cou=co-up%*%x
if(is.null(mnorm)) {
mban=ban
if (is.null(x0)) {
x0=-x
} else {
x0=x0-x
}
} else {
mban=mnorm%*%ban
if (is.null(x0)) {
x0=-mnorm%*%x
} else {
x0=mnorm%*%(x0-x)
}
}
z=lsi_ln(mban, x0, uz, cou)
if (!anyNA(z)) {
x=x+ban%*%z
attr(x, "mes")=join("\n", attr(x, "mes"), attr(z, "mes"))
} else {
#print(list(a=a, b=b, u=u, co=co))
#browser()
attr(x, "mes")=join("\n", attr(x, "mes"), "lsi_ln: Unfeasible constraints in null space. (It must not be. Round off errors struck again.)")
}
}
}
x[aq$pivot]=x
if (!is.null(u)) {
cou=co-u%*%x
if (any(cou > tol)) {
# round errors made fail inequality constraints
# force them as much as possible even by degrading the residual
dx=ldp(u, cou, rcond)
if (!is.null(dx)) {
x=x+dx
} # else leave x as is
}
#stopifnot(all(u%*%x-co >= -tol))
}
attr(x, "mes")=join("\n", attr(x, "mes"), paste(
"lsi_ln: Rank deficient matrix in least squares\n",
ndefic, " free variable(s):\n",
paste(dimnames(a)[[2]][aq$pivot[i2]],
aq$pivot[-i1], sep="\t", collapse="\n"),
"\nLeast L2-norm solution is provided.",
sep=""))
return(x)
}
#' Least Distance Problem
#'
#' @description
#' Solve least distance programming: find x satisfying \code{u\%*\%x >= co} and s.t. min(||x||)
#' by passing to nnls (non negative least square) problem.
#' @param u a dense matrix of inequality constraints
#' @param co a right hand side vector of inequality constraints
#' @param rcond maximal condition number for determining rank deficient matrix
#' @return solution vector or NULL if constraints are unfeasible
#' @importFrom nnls nnls
#' @export
ldp=function(u, co, rcond=1.e10) {
m=NROW(u)
n=ncol(u)
tol=1.e-10
if (m == 0) {
# no inequality to satisfy => trivial case
return(double(n))
}
rcond=abs(rcond)
# eliminate NA from co
i0=which(is.na(co))
if (length(i0) > 0) {
u=u[-i0,,drop=FALSE]
co=co[-i0]
m=m-length(i0)
if (m == 0) {
# no inequality to satisfy => all NA
return(rep(NA_real_, n))
}
}
# eliminate 0 rows from u
maxu=apply(u, 1, function(v) max(abs(v)))
i0=which(maxu < 1.e-13)
if (length(i0) > 0) {
if (all(co[i0] < tol)) {
u=u[-i0,,drop=FALSE]
co=co[-i0]
m=m-length(i0)
if (m == 0) {
# no inequality to satisfy => trivial case
return(double(n))
}
} else {
# unfeasible 0 constraints
return(NULL)
}
}
maxco=max(co)
if (maxco<=tol) {
# all rhs are <= 0 => trivial case
return(double(n))
}
e=rbind(t(u), t(co))
f=double(n+1L)
f[n+1]=1.
resnnls=nnls::nnls(e, f)
feasible=sqrt(resnnls$deviance) > tol && resnnls$residuals[n+1] != 0.
if (feasible) {
x=resnnls$residuals[1:n]/(-resnnls$residuals[n+1])
# check for numerical stability problems
ux=u%*%x
cou=co-ux
if (FALSE && any(cou > tol)) { # FALSE because it can worsen the situation
# second trial
e[n+1,]=cou
rn=nnls::nnls(e, f)
if (rn$residuals[n+1]!=0.) {
x=x-rn$residuals[1:n]/rn$residuals[n+1]
} else {
x=NULL
}
} else if (all(cou<0) && !all(x==0.)) {
# round off error pushed the solution inside of feasible domain
# shorten it till the closest border
alpha=ux/co
alpha=min(alpha[alpha>=1.], na.rm=T)
x=x/alpha
}
#stopifnot(all(u%*%x-co >= -tol))
} else {
x=NULL
}
return(x)
}
norm2=function(v)crossprod(as.numeric(v))[1L]
ls_ln_old=function(a, b, rcond=1.e10) {
# least squares with least norm
# no LAPACK in the second qr()
if (! is.qr(a)) {
a=base::qr(a, LAPACK=TRUE)
}
n=ncol(a$qr)
d=diag(a$qr)
a$rank=sum(abs(d)>abs(d[1])/rcond)
rdefic=a$rank < n
i1=1:a$rank
i2=(a$rank+1):n
r=base::qr.R(a)[i1,,drop=FALSE]
if (!rdefic) {
# plain ls
x=base::backsolve(r, qr.qty(a, b))
x[a$pivot]=x
return(x)
}
# prepare free variable substitution
r1=r[,i1,drop=FALSE]
r2=r[,i2,drop=FALSE]
qrr=qr(t(r), LAPACK=F)
# null space and its complement bases
ba=base::qr.Q(qrr, complete=TRUE)
bal=ba[,i1,drop=FALSE]
# least norm
x=bal%*%base::backsolve(qr.R(qrr), qr.qty(a, b), transpose=TRUE)
x[a$pivot]=x
return(x)
}
#' Linear Least Squares, least norm solution
#'
#' @param a matrix or its QR decomposition
#' @param b vector of right hand side
#' @param rcond maximal condition number for rank definition
#' @return solution vector
#' @export
ls_ln=function(a, b, rcond=1.e10) {
# least squares with least norm
# LAPACK is called in the second qr() too.
if (! is.qr(a)) {
a=base::qr(a, LAPACK=TRUE)
}
n=ncol(a$qr)
d=diag(a$qr)
a$rank=sum(abs(d)>abs(d[1])/rcond)
if (a$rank == 0) {
return(double(n)) # return plain 0 vector
}
rdefic=a$rank < n
i1=seq(len=a$rank)
r=base::qr.R(a)[i1,,drop=FALSE]
if (!rdefic) {
# plain ls
x=base::backsolve(r, qr.qty(a, b))
x[a$pivot]=x
return(x)
}
ndefic=n-a$rank
i2=seq(a$rank+1, n, len=ndefic)
# prepare free variable substitution
qrr=base::qr(t(r), LAPACK=TRUE)
# least norm
x=base::qr.qy(qrr, c(backsolve(qrr$qr, qr.qty(a, b)[qrr$pivot], transpose=TRUE), rep.int(0., ndefic)))
x[a$pivot]=x
return(x)
}
#' Linear Least Squares problem with inequality and equality constraints, least norm solution
#'
#' @description
#' solve linear least square problem (min ||A%*%x-b||)
#' with inequality constraints u%*%x>=co and equality constraints e%*%x=ce
#' Method:
#' reduce the pb to lsi_ln on the null-space of e
#' @param a dense matrix A or its QR decomposition
#' @param b right hand side vector
#' @param u dense matrix of inequality constraints
#' @param co right hand side vector of inequality constraints
#' @param e dense matrix of equality constraints
#' @param ce right hand side vector of equality constraints
#' @param rcond maximal condition number for determining rank deficient matrix
#' @return solution vector whose attribute 'mes' may contain a message about possible numerical problems
#' @seealso [lsi_ln]
#' @export
lsie_ln=function(a, b, u=NULL, co=NULL, e=NULL, ce=NULL, rcond=1.e10) {
if (is.qr(a)) {
n=ncol(a$qr)
aqr=TRUE
} else {
n=ncol(a)
aqr=FALSE
}
if (!is.null(e) && nrow(e) > 0) {
# x is searched in a form xp+bn*z
# where xp is a particular solution of e*x=ce
# bn is a basis of e null space and z is a
# new unknown vector.
nr=nrow(e)
nc=ncol(e)
if (nr > nc) {
x=rep(NA, n)
attr(x, "mes")="lsie_ln: matrix e is over determined."
return(x)
}
bn=Nulla(t(e))
qe=attr(bn, "qr")
if (qe$rank < nr) {
x=rep(NA, n)
attr(x, "mes")="lsie_ln: matrix e is rank deficient."
return(x)
}
# particular solution
xp=base::qr.qy(qe, c(backsolve(qe$qr, ce[qe$pivot], transpose=T), double(n-nr)))
if (aqr) {
R=base::qr.R(a, complete=TRUE)
b=b-base::qr.qy(a, R%*%xp[a$pivot])
a=base::qr.qy(a, R%*%bn[a$pivot,,drop=FALSE])
} else {
b=b-a%*%xp
a=a%*%bn
}
if (!is.null(u)) {
co=co-u%*%xp
u=u%*%bn
}
z=lsi_ln(a, b, u, co, rcond)
x=xp+bn%*%z
if (!is.null(attr(z, "mes"))) {
attr(x, "mes")=attr(z, "mes")
}
x
} else {
lsi_ln(a, b, u=u, co=co, rcond=rcond)
}
}
#' Linear Least Squares, least norm solution (by svd)
#'
#' @description
#' Least squares \code{a\%*\%x ~= b} of least norm ||x|| by using svd(a)
#' @param a dense matrix
#' @param b right hand side vector
#' @param rcond maximal condition number for determining rank deficient matrix
#' @return solution vector
#' @export
ls_ln_svd=function(a, b, rcond=1.e10) {
sa=svd(a)
d=sa$d
rank=sum(d > d[1L]/rcond)
if (rank == 0L) {
return(rep(0., ncol(a)))
}
i=seq_len(rank)
x=sa$v[,i]%*%(crossprod(sa$u[,i], b)/d[i])
return(x)
}
#' Total Least Squares \code{a\%*\%x ~= b}
#'
#' @param a matrix
#' @param b right hand side vector
#' @return solution vector
#' @export
tls=function(a, b) {
# total least square by svd
sab=svd(cbind(a, b))
n=ncol(a)
v=sab$v[,length(sab$d)]
return(v[seq_len(n)]/(-v[n+1L]))
}
#' Regularized Linear Least Squares
#'
#' @description
#' solve linear least square problem \code{(min_x ||a*x-b||)}
#' with inequality constraints ux>=co
#' If a is rank deficient, regularization term \code{lambda^2*||mnorm*(x-x0)||^2}
#' is added to \code{||a*x-b||^2}.
#' @details
#' The rank of a is estimated as number of singular values
#' above \code{d[1]*1.e-10} where \code{d[1]} is the highest singular value.
#' The scalar lambda is an positive number
#' and is calculated as \code{d[1]/sqrt(rcond)} ('rcond' parameter is preserved
#' for compatibility with others lsi_...() functions).
#' At return, lambda can be found in attributes of the returned vector x.
#' NB. lambda is set to NA
#' - if rank(a)==0 or a is of full rank
#' - or if there is no inequality.
#' If the matrix mnorm is NULL, it is supposed to be an identity matrix.
#' If the vector x0 is NULL, it is treated as 0 vector.
#' @param a dense matrix A or its QR decomposition
#' @param b right hand side vector
#' @param u dense matrix of inequality constraints
#' @param co right hand side vector of inequality constraints
#' @param rcond used for calculating \code{lambda=d[1]/sqrt(rcond)} where \code{d[1]} is maximal singular value of a
#' @param mnorm norm matrix (can be dense or sparse) for which %*% operation with a dense vector is defined
#' @param x0 optional vector from which a least norm distance is searched for
#' @return solution vector whose attribute 'mes' may contain a message about possible
#' numerical problems and 'lambda' is regularization parameter used in solution.
#' @seealso [lsi_ln]
#' @export
lsi_reg=function(a, b, u=NULL, co=NULL, rcond=1.e10, mnorm=NULL, x0=NULL) {
mes=""
lambda=1./sqrt(rcond)
tol=1.e-10
if (is.qr(a)) {
# get back the matrix from qr
a=base::qr.X(a)
}
n=ncol(a)
m=nrow(a)
if (m < n) {
# make nrow(a) >= ncol(a) by repeating a and b as many times as needed
nrep=ceiling(n/m)
a=aperm(array(a, dim=c(m, n, nrep)), c(1,3,2))
dim(a)=c(m*nrep, n)
b=rep(b, nrep)
}
sva=svd(a)
d=sva$d
arank=sum(d > d[1]*tol)
x0=if (is.null(x0)) double(n) else x0
if (arank == 0) {
# matrix a is zero => no x can diminish the residual norm
if (!is.null(u) && nrow(u) > 0) {
# there are some inequalities to deal with
co0=co-u%*%x0
if (any(co0 > tol)) {
# at least we satisfy inequalities
if (is.null(mnorm)) {
dx=ldp(u, co0)
if (is.null(dx)) {
x=rep(NA, n)
attr(x, "mes")="lsi_reg: matrix a is zero rank and ldp() revealed unfeasible constrants"
}
} else {
x=lsi_reg(mnorm, mnorm%*%x0, u, co, rcond, NULL, NULL)
}
} else {
# inequalities are already satisfied
x=x0
}
} else {
# no inequalities
x=x0
}
attr(x, "lambda")=NA
return(x)
}
# hereafter, the arank is > 0
i=seq_len(arank)
if (NROW(u) == 0) {
# no inequality
# plain ls_ln
if (arank == n) {
x=sva$v[,i]%*%(crossprod(sva$u[,i,drop=FALSE], b)/d[i])
attr(x, "lambda")=NA
} else {
mes=sprintf("lsi_reg: Rank deficient matrix in least squares
There is (are) %d free variable(s).
Regularized L2-norm solution is provided.", n-arank)
if (is.null(mnorm)) {
mnorm=diag(1., n)
}
lambda=d[1]*lambda
lam2=lambda**2
d2=diag(d**2, n)
mv=mnorm%*%sva$v
mv2=crossprod(mv)
bx0=crossprod(mv, mnorm%*%x0)
bt=d*crossprod(sva$u, b)+lam2*x0
x=sva$v%*%solve(d2+lam2*mv2, d*bt)
attr(x, "lambda")=lambda
attr(x, "mes")=mes
}
} else {
# there are inequalities
if (arank == n) {
# a is of full rank => no lambda to use, plain lsi
invd=1./d
bt=crossprod(sva$u, b)
ut=(u%*%sva$v)*rep(invd, rep(n, n))
y=ldp(ut, co-ut%*%bt)
x=sva$v%*%(invd*(y+bt))
attr(x, "lambda")=NA
} else {
mes=sprintf("lsi_reg: Rank deficient matrix in least squares
There is (are) %d free variable(s).
Regularized L2-norm solution is provided.", n-arank)
if (is.null(mnorm)) {
mnorm=diag(1., n)
}
# transformed rhs
mv=mnorm%*%sva$v
bx0=crossprod(mv, mnorm%*%x0)
bt=d*crossprod(sva$u, b)
# semi-transformed u (it lacks a factor 1/(d**2+(lambda*mv)**2) )
uv=u%*%sva$v
mv2=crossprod(mv)
# ldp to find approximate solution
lambda=d[1]*lambda
lam2=lambda**2
d2=diag(d**2, n)+lam2*mv2
ut=t(solve(t(d2), t(uv)))
btt=bt+lam2*bx0
y=ldp(ut, co-ut%*%btt)
x=sva$v%*%solve(d2, btt+y)
attr(x, "lambda")=lambda
attr(x, "mes")=mes
}
}
return(x)
}
#' Transform box-type inequalities into matrix and vector form
#'
#' @description
#' Transform a set of inequalities \code{param["name"] >= lower["name"]} and
#' \code{param["name"] <= upper["name"]} into a list with matrix u and vector co such that
#' u%*%param>=co. In addition to box inequalities above, user can provide linear
#' inequalities in a form like \code{"a+2*c+3*b >= 0"} where 'a', 'b' and 'c' must be names of
#' param components. Numeric and symbolic coefficients and right hand sides
#' are allowed in these expressions. However, symbols must be defined at the moment
#' of calling \code{uplo2uco()} so that expressions containing such symbols could
#' be \code{eval()}-ed to numerical values. All inequalities must be written with '>='
#' sign (not with '<=', '>', ...).
#' @param param a named vector whose names will be used for parsing inequalities
#' @param upper a named numeric vector of upper limits
#' @param lower a named numeric vector of lower limits
#' @param linear a string vector of linear inequalities
#' @return a list with numeric matrix 'u' and vector 'co' such that \code{u%*%param-co>=0}
#' @seealso [equa2vecmat] for parsing linear expressions
#' @export
uplo2uco=function(param, upper=NULL, lower=NULL, linear=NULL) {
nm_par=names(param)
nlo=length(lower)
nup=length(upper)
nli=length(linear)
u=matrix(0., nrow=nup+nlo+nli, ncol=length(param))
co=numeric(nrow(u))
colnames(u)=nm_par
rownames(u)=c(
if (nup) paste(names(upper), " <= ", upper, sep="") else NULL,
if (nlo) paste(names(lower), " >= ", lower, sep="") else NULL,
linear
)
names(co)=rownames(u)
# fill u and co
if (nup) {
u[seq_len(nup), names(upper)]=diag(-1, nup)
co[seq_len(nup)]=-upper
}
u[nup+seq_len(nlo), names(lower)]=diag(1, nlo)
co[nup+seq_len(nlo)]=lower
# parse inequalities
if (nli) {
vema=try(equa2vecmat(nm_par, linear, sep=">="))
if (inherits(vema, "try-error")) {
stop("Error in parsing inequalities")
}
i=nup+nlo+seq_len(nli)
u[i, ]=vema[,-1L]
co[i]=vema[,1L]
}
return(list(u=u, co=co))
}
#' Parse linear equations/inequalities
#'
#' @description
#' parse a text vector of linear equations and produce a corresponding
#' matrix and right hand side vector
#' @param nm_par a string vector of variable names. It will be used in the symbolic
#' derivation.
#' @param linear string vector of linear equations like \code{"a+2*c+3*b = 0"}
#' @param sep separator of two parts of equations. Use for example
#' ">=" for linear inequalities
#' @return an augmented matrix. Its first column is the rhs vector.
#' Other columns are named by nm_par. If the vector linear is NULL or its content
#' is empty a NULL is returned
#' @examples
#' equa2vecmat(c("a", "b", "c"), "a+2*c+3*b = 0", "=")
#' @export
equa2vecmat=function(nm_par, linear, sep="=") {
# Sanity check
stopifnot(length(sep)==1 && nchar(sep) > 0)
stopifnot(length(nm_par)>=1 && all(nchar(nm_par)) > 0)
vlin=sapply(linear, function(it) strsplit(it, sep)[[1]])
if (length(vlin) && !is.matrix(vlin)) {
stop(sprintf("Linear (in)equalities are expected to have '%s' sign in them", sep))
}
if (length(vlin) == 0) {
# we are set
return(NULL)
}
# each column in vlin is an (in)equality
# first row is left part, the second row is the right part
# derive left and right parts to get the matrix
ze=double(length(nm_par))
names(ze)=nm_par
ze=as.list(ze)
de=apply(vlin, 2L, function(ineq) {
le=with(ze, eval(stats::deriv(parse(text=ineq[1]), nm_par)))
ri=with(ze, eval(stats::deriv(parse(text=ineq[2]), nm_par)))
v=c(ri-le, attr(le, "gradient")-attr(ri, "gradient"))
return(v)
})
rownames(de)=c("rhs", nm_par)
return(t(de))
}
lsi_lim=function(a, b, u=NULL, co=NULL, rcond=1.e10, mnorm=NULL, x0=NULL) {
if (requireNamespace("limSolve", quietly = TRUE))
suppressWarnings(limSolve::lsei(A=a, B=b, G=u, H=co)$X)
}
#' Null-space basis
#'
#' @description
#' use Lapack for null space basis
#' (derived from MASS::Null)
#' @param M matrix such that \code{t(M)\%*\%B=0} where B is a basis of t(M)'s kernel (aka Null-space)
#' @param rcond maximal condition number for rank definition
#' @return numeric matrix whose columns are basis vectors. Its attribute 'qr' contains QR decomposition of M.
#' @examples
#' Nulla(1:3)
#' @seealso [MASS::Null]
#' @export
Nulla=function (M, rcond=1.e10) {
tmp <- qr(as.matrix(M), LAPACK=TRUE)
d=abs(diag(tmp$qr))
n=length(d)
#browser()
if (d[1L]==0.) {
tmp$rank = 0
} else {
tmp$rank = sum(d/d[1L] > 1./rcond)
}
set <- if (tmp$rank == 0L)
1L:n
else
-(1L:tmp$rank)
bn=qr.Q(tmp, complete = TRUE)[, set, drop = FALSE]
attr(bn, "qr")=tmp
return(bn)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.