| Title: | Non Linear Least Squares with Inequality Constraints |
|---|---|
| Description: | We solve non linear least squares problems with optional equality and/or inequality constraints. Non linear iterations are globalized with back-tracking method. Linear problems are solved by dense QR decomposition from 'LAPACK' which can limit the size of treated problems. On the other side, we avoid condition number degradation which happens in classical quadratic programming approach. Inequality constraints treatment on each non linear iteration is based on 'NNLS' method (by Lawson and Hanson). We provide an original function 'lsi_ln' for solving linear least squares problem with inequality constraints in least norm sens. Thus if Jacobian of the problem is rank deficient a solution still can be provided. However, truncation errors are probable in this case. Equality constraints are treated by using a basis of Null-space. User defined function calculating residuals must return a list having residual vector (not their squared sum) and Jacobian. If Jacobian is not in the returned list, package 'numDeriv' is used to calculated finite difference version of Jacobian. The 'NLSIC' method was fist published in Sokol et al. (2012) <doi:10.1093/bioinformatics/btr716>. |
| Authors: | Serguei Sokol [aut, cre] |
| Maintainer: | Serguei Sokol <[email protected]> |
| License: | GPL-2 |
| Version: | 1.0.4 |
| Built: | 2025-02-25 04:45:21 UTC |
| Source: | https://github.com/mathscell/nlsic |
parse a text vector of linear equations and produce a corresponding matrix and right hand side vector
equa2vecmat(nm_par, linear, sep = "=")equa2vecmat(nm_par, linear, sep = "=")
nm_par |
a string vector of variable names. It will be used in the symbolic derivation. |
linear |
string vector of linear equations like |
sep |
separator of two parts of equations. Use for example ">=" for linear inequalities |
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
equa2vecmat(c("a", "b", "c"), "a+2*c+3*b = 0", "=")equa2vecmat(c("a", "b", "c"), "a+2*c+3*b = 0", "=")
convert elements of vector v (and all following arguments) in strings and join them using sep as separator.
join(sep, v, ...)join(sep, v, ...)
sep |
A string used as a separator |
v |
A string vector to be joined |
... |
other variables to be converted to strings and joined |
A joined string
join(" ", c("Hello", "World"))join(" ", c("Hello", "World"))
Solve least distance programming: find x satisfying u%*%x >= co and s.t. min(||x||)
by passing to nnls (non negative least square) problem.
ldp(u, co, rcond = 1e+10)ldp(u, co, rcond = 1e+10)
u |
a dense matrix of inequality constraints |
co |
a right hand side vector of inequality constraints |
rcond |
maximal condition number for determining rank deficient matrix |
solution vector or NULL if constraints are unfeasible
Linear Least Squares, least norm solution
ls_ln(a, b, rcond = 1e+10)ls_ln(a, b, rcond = 1e+10)
a |
matrix or its QR decomposition |
b |
vector of right hand side |
rcond |
maximal condition number for rank definition |
solution vector
Least squares a%*%x ~= b of least norm ||x|| by using svd(a)
ls_ln_svd(a, b, rcond = 1e+10)ls_ln_svd(a, b, rcond = 1e+10)
a |
dense matrix |
b |
right hand side vector |
rcond |
maximal condition number for determining rank deficient matrix |
solution vector
solve linear least square problem (min ||A%%x-b||)
with inequality constraints u%%x>=co
lsi(a, b, u = NULL, co = NULL, rcond = 1e+10, mnorm = NULL, x0 = NULL)lsi(a, b, u = NULL, co = NULL, rcond = 1e+10, mnorm = NULL, x0 = NULL)
a |
dense matrix A or its QR decomposition |
b |
right hand side vector. Rows containing NA are dropped. |
u |
dense matrix of inequality constraints |
co |
right hand side vector of inequality constraints |
rcond |
maximal condition number for determining rank deficient matrix |
mnorm |
dummy parameter |
x0 |
dummy parameter |
Method:
reduce the problem to ldp (min(xat*xa) => least distance programming)
solve ldp
change back to x
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
solution vector whose attribute 'mes' may contain a message about possible numerical problems
solve linear least square problem min_x ||A*x-b||
with inequality constraints u%*%x >= co
If A is rank deficient, least norm solution ||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.
lsi_ln(a, b, u = NULL, co = NULL, rcond = 1e+10, mnorm = NULL, x0 = NULL)lsi_ln(a, b, u = NULL, co = NULL, rcond = 1e+10, mnorm = NULL, x0 = NULL)
a |
dense matrix A or its QR decomposition |
b |
right hand side vector |
u |
dense matrix of inequality constraints |
co |
right hand side vector of inequality constraints |
rcond |
maximal condition number for determining rank deficient matrix |
mnorm |
norm matrix (can be dense or sparse) for which |
x0 |
optional vector from which a least norm distance is searched for |
solution vector whose attribute 'mes' may contain a message about possible numerical problems
solve linear least square problem (min_x ||a*x-b||)
with inequality constraints ux>=co
If a is rank deficient, regularization term lambda^2*||mnorm*(x-x0)||^2
is added to ||a*x-b||^2.
lsi_reg(a, b, u = NULL, co = NULL, rcond = 1e+10, mnorm = NULL, x0 = NULL)lsi_reg(a, b, u = NULL, co = NULL, rcond = 1e+10, mnorm = NULL, x0 = NULL)
a |
dense matrix A or its QR decomposition |
b |
right hand side vector |
u |
dense matrix of inequality constraints |
co |
right hand side vector of inequality constraints |
rcond |
used for calculating |
mnorm |
norm matrix (can be dense or sparse) for which %*% operation with a dense vector is defined |
x0 |
optional vector from which a least norm distance is searched for |
The rank of a is estimated as number of singular values
above d[1]*1.e-10 where d[1] is the highest singular value.
The scalar lambda is an positive number
and is calculated as 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.
solution vector whose attribute 'mes' may contain a message about possible numerical problems and 'lambda' is regularization parameter used in solution.
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
lsie_ln(a, b, u = NULL, co = NULL, e = NULL, ce = NULL, rcond = 1e+10)lsie_ln(a, b, u = NULL, co = NULL, e = NULL, ce = NULL, rcond = 1e+10)
a |
dense matrix A or its QR decomposition |
b |
right hand side vector |
u |
dense matrix of inequality constraints |
co |
right hand side vector of inequality constraints |
e |
dense matrix of equality constraints |
ce |
right hand side vector of equality constraints |
rcond |
maximal condition number for determining rank deficient matrix |
solution vector whose attribute 'mes' may contain a message about possible numerical problems
Solve non linear least squares problem min_par ||r(par,...)$res||
with optional inequality constraints u%*%par >= co
and optional equality constraints e%*%par = eco
nlsic( par, r, u = NULL, co = NULL, control = list(), e = NULL, eco = NULL, flsi = lsi, ... )nlsic( par, r, u = NULL, co = NULL, control = list(), e = NULL, eco = NULL, flsi = lsi, ... )
par |
initial values for parameter vector. It can be in non feasible domain,
i.e. for which |
r |
a function calculating residual vector
by a call to |
u |
constraint matrix in |
co |
constraint right hand side vector |
control |
a list of control parameters ('=' indicates default values):
|
e |
linear equality constraint matrix in |
eco |
right hand side vector in equality constraints |
flsi |
function solving linear least squares problem. For a custom function,
see interfaces in |
... |
parameters passed through to r() |
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 lsi() call.
NB. If the function 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.
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().
NB3. User should test that field 'mes' is not NULL even when error is 0. It may
contain a warning message.
a list with following components (some components can be absent depending on 'control' parameter)
'par' estimated values of par
'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 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
# 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) } }# 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) } }
use Lapack for null space basis (derived from MASS::Null)
Nulla(M, rcond = 1e+10)Nulla(M, rcond = 1e+10)
M |
matrix such that |
rcond |
maximal condition number for rank definition |
numeric matrix whose columns are basis vectors. Its attribute 'qr' contains QR decomposition of M.
Nulla(1:3)Nulla(1:3)
a%*%x ~= b
Total Least Squares a%*%x ~= b
tls(a, b)tls(a, b)
a |
matrix |
b |
right hand side vector |
solution vector
Transform a set of inequalities param["name"] >= lower["name"] and
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 "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 uplo2uco() so that expressions containing such symbols could
be eval()-ed to numerical values. All inequalities must be written with '>='
sign (not with '<=', '>', ...).
uplo2uco(param, upper = NULL, lower = NULL, linear = NULL)uplo2uco(param, upper = NULL, lower = NULL, linear = NULL)
param |
a named vector whose names will be used for parsing inequalities |
upper |
a named numeric vector of upper limits |
lower |
a named numeric vector of lower limits |
linear |
a string vector of linear inequalities |
a list with numeric matrix 'u' and vector 'co' such that u%*%param-co>=0
equa2vecmat for parsing linear expressions