LP within get_lambda_max (silently) fails converge to a solution

[huisaddison]

Reproducing example:

library(quantgen)

set.seed(1)
n = 100
p = 30
d = get_diff_mat(p, 3)

X = matrix(rnorm(n*p), nrow=n)
y = rnorm(n)

get_lambda_max(X, y, d=d) # Returns 72.50772

get_lambda_max(X, y, d=d, lp_solver='gurobi') # Returns NULL

If one inspects the returned solution object a following the lines

quantgen/quantgen/R/quantile_genlasso.R

Lines 557 to 573 in 5971854

	# Gurobi
	if (use_gurobi) {
	model$sense = c(rep(">=", 2*m), rep("=", p))
	model$lb = c(rep(-Inf,m), 0)
	a = gurobi::gurobi(model=model, params=list(LogToConsole=0))
	lambda_max = a$x[m+1]
	}
	# GLPK
	else {
	model$sense = c(rep(">=", 2*m), rep("==", p))
	model$bounds = list(lower=list(ind=1:m, val=rep(-Inf,m)))
	a = Rglpk_solve_LP(obj=model$obj, mat=model$A, dir=model$sense,
	rhs=model$rhs, bounds=model$bounds)
	lambda_max = a$solution[m+1]
	}

one will find that with GLPK,

Called from: get_lambda_max(X, y, d = d)
Browse[1]> print(a$status)
[1] 1

and with Gurobi,

Called from: get_lambda_max(X, y, d = d, lp_solver = "gurobi")
Browse[1]> print(a$status)
[1] "INFEASIBLE"

The failure is more obvious under Gurobi because the solution will be NULL, whereas with Rglpk the supposed solution is the final iterate before it quits.

Thanks @jeremy-goldwasser for catching this

[huisaddison]

Next steps:

• Is the implementation correct? (i.e., are we forming the LP determining lambda_max correctly?)
• And if it is correct, do we expect Gurobi/GLPK to fail to find the solution often? What safety logic should we put in / heuristics can we use instead for a default large value of lambda?

[huisaddison]

Just adding some notes to this so we don't forget (last discussed on 22-Dec-2022):

• Current implementation is based on a heuristic that we don't have documentation for; and was only ever meant to give an approximate value of lambda_max.
• We can replace it with a different kind of heuristic, e.g., start with lambda_max in the least squares case (which can be obtained in closed form), apply an inflation factor to account for the least absolute deviations case. And maybe "explore", doubling lambda_max until we find that the active set is empty (check the empty active set against the KKT conditions)