# Nonnegative identity-link Poisson regression with ridge or fused ridge penalty

I would like to fit nonnegative identity-link Poisson regression models with a ridge or fused ridge penalty, i.e. with nonnegativity constraints on the fitted coefficients, Poisson error noise & a ridge or fused ridge penalty on the fitted coefficients.

In the case of standard nonnegative least squares ridge regression I understand that this can be fit very easily by row augmenting the original covariate matrix X with a diagonal p x p matrix with sqrt(lambdas) along the diagonal and augmenting the outcome vector y with p zeros and then doing a nonnegative least squares (nnls) regression using these augmented matrices/vectors.

Two question though:

1) Does this recipe also work for nonnegative identity-link Poisson models, the idea being that I would then simply change from an nnls fit to a nnpois fit (fit using an EM algorithm, using the addreg package). Would anybody be able to prove analytically that this would indeed be equivalent to optimizing the ridge penalized nonnegative identity-link Poisson likelihood?

2) Would a recipe like this also work to fit nonnegative least square or nonnegative Poisson models with a fused ridge penalty, where the row augmented part of the covariate matrix would then no longer be diag(sqrt(lambda),p)) as in regular ridge regression, but sqrt(lambda)*(t(D) %*% D) with D=diff(diag(p), 1). If I try this, this works if I don't impose nonnegativity constraints (e.g. with a regular lm.fit), but it fails as soon as I plug this into either nnls or nnpois, as they only take covariate matrices with all positive elements. Is there an easy fix for this?

Regarding question 1) I did some quick numerical tests and I do indeed get a very similar results based on using nnpois with a covariate matrix row augmented with diag(sqrt(lambda),p)) and based on directly optimizing the penalized ML objective, but I am on the lookout for an analytical proof that this is a correct thing to do (or if it is not exact, that it would at least be a very close approximation of the actual penalized ML objective)... Any thoughts?

Numerical example, here in the context of trying to estimate a peak shape function from a superposition of blurred impulses with known pulse locations but unknown blur kernel, and where the noise on the measured superposition of the blurred peaks would be Poisson:

# simulated data
require(Matrix)
n = 200
x = 1:n
npeaks = 20
set.seed(123)
u = sample(x, npeaks, replace=FALSE) # peak locations, which are assumed to be known here
peakhrange = c(10,1E3) # peak height range
h = 10^runif(npeaks, min=log10(min(peakhrange)), max=log10(max(peakhrange))) # peak heights, which are assumed to be known here
a = rep(0, n) # locations of spikes of simulated spike train, which are assumed to be known here
a[u] = h
gauspeak = function(x, u, w, h=1) h*exp(((x-u)^2)/(-2*(w^2))) # simulated peak shape, assumed to be unknown
bM = do.call(cbind, lapply(1:n, function (u) gauspeak(x, u=u, w=5, h=1) )) # banded matrix with theoretical peak shape function used
p = 50 # desired size of peak shape model
g_theor = gauspeak(x=1:p, u=p/2, w=5, h=1) # theoretical unknown peak shape function
y_nonoise = as.vector(bM %*% a) # noiseless simulated signal = linear convolution of spike train with peak shape function
y = rpois(n, y_nonoise) # simulated signal with random poisson noise on it
G = matrix(0, n + p - 1, p)
G = Matrix(G, sparse=TRUE)
for (k in 1:p) G[(1:n) + k - 1, k] = a
G = G[(1:n) + floor(p/2) - 1, ]
par(mfrow=c(1,1))
plot(y, type="l", ylab="Signal", xlab="x", main="Known spike train (red) convoluted with unknown blur kernel and Poisson noise")
lines(a, type="h", col="red")


# peak shape function estimated using weighted nnls ridge regression via nnls with a row augmented covariate matrix
weights = 1/(y+1) # 1/variance obs weights to account for Poisson noise in weighted nonnegative least square approximation
weights = n*weights/sum(weights)
library(nnls)
lambda = 1E4 # regularization parameter for ridge penalty
library(microbenchmark)
microbenchmark(g_wnnls <- nnls(A=as.matrix(rbind(G*sqrt(weights), diag(sqrt(lambda),p))), b=c(y*sqrt(weights),rep(0,p)))$x) # 1.4 ms g_wnnls = (g_wnnls-min(g_wnnls))/(max(g_wnnls)-min(g_wnnls)) dev.off() par(mfrow=c(2,1)) plot(g_theor, type="l", lwd=7, col="grey", main="Ground truth (grey), nonnegative weighted LS ridge estimate (red)", ylab="Peak shape", xlab="x") lines(g_wnnls, col="red", lwd=2) # peak shape function estimated using nonnegative Poisson adaptive ridge regression via nnpois with a row augmented covariate matrix library(addreg) lambda = 1E6 # regularization parameter for ridge penalty adpenweights = (1/(g_wnnls^2+1E-5)) # adaptive penalty weights to do adaptive ridge regression using wnnls estimates as truncated Gaussian prior adpenweights = n*adpenweights/sum(adpenweights) lambdas = lambda*adpenweights # adaptive penalties for adaptive ridge penalty system.time(g_nnpoisridge <- nnpois(y=c(y,rep(0,p)), # nonnegative identity link Poisson adaptive ridge regression, solved using EM algo x=as.matrix(rbind(G, diag(sqrt(lambdas),p))), standard=rep(1,p), offset=0, start=rep(1E-5, p), # or we can use g_wnnls+1E-5 accelerate="squarem")$coefficients) # 0.02
g_nnpoisridge = (g_nnpoisridge-min(g_nnpoisridge))/(max(g_nnpoisridge)-min(g_nnpoisridge))
plot(gauspeak(x=1:p, u=25, w=5, h=1), type="l", lwd=7, col="grey", main="Ground truth (grey), nonnegative Poisson adaptive ridge estimate (red)", ylab="Peak shape", xlab="x")
lines(g_nnpoisridge, col="red", lwd=2)
# PS: instead of using nnpois we could also use nnlm (NNLM package) with loss="mkl" (Kullback-Leibler) which is then solved using coordinate descent and gives almost the same result


# check that nonnegative poisson adaptive ridge regression is calculated correctly using ML maximization with port algo (Quasi-Newton BFGS)
# negative log-likelihood for identity link poisson ridge regression
NLL_RIDGEPOIS = function(coefs, X, y, lambdas) {
preds = X %*% coefs
LLs = stats::dpois(y, lambda=preds, log=TRUE) # log likelihood contributions of each observation
n = nrow(X)
-(1/n)*sum(LLs) - lambdas*sum(coefs^2)
}
g_nnpoisridgeML = nlminb(start = g_nnpoisridge, # we use our estimates above as starting values
objective = NLL_RIDGEPOIS, X = as.matrix(G), y = y, lambdas=lambdas,
control = list(iter.max=1000, abs.tol=1E-20, rel.tol=1E-15),
lower = rep(0,p) # we use nonnegativity constraints
)\$par
max(g_nnpoisridgeML) # not sure why I get such a large value here, 6.269765e+68, should be around 1 as above - maybe mistake here?
# as a quick fix I rescale coefficients
g_nnpoisridgeML = (g_nnpoisridgeML-min(g_nnpoisridgeML))/(max(g_nnpoisridgeML)-min(g_nnpoisridgeML))
max(abs(g_nnpoisridgeML-g_nnpoisridge)) # 0.01726672 - ML estimate is close to solution above, but not identical - could be optimizer issue though or mistake in code...
par(mfrow=c(1,1))
plot(g_theor,col="grey", type="l", lwd=5, main="Ground truth (grey), nonnegative Poisson adaptive ridge estimate (red),\n nonnegative Poisson adaptive ridge ML estimate (blue)", ylab="Peak shape", xlab="x")
lines(g_nnpoisridge,col="red",lwd=2)
lines(g_nnpoisridgeML,col="blue",lwd=2)


EDIT: Just figured that in the end it may be exact - reason is that you can readily fit nonnegative GLMs using iteratively reweighted nonnegative least squares (just like the regular IRLS algo but substituting the regular solve / lm.wfit by weighted nnls) and for identity link you just have to regress yhat on your covariate matrix and the weights W are just 1/expected variance (ie in the IRLS algo, http://bwlewis.github.io/GLM/, z=yhat and W=1/expected variance=1/yhat for Poisson noise).

So given that we are just doing weighted (nonnegative) least squares anyway, adding a ridge penalty in each iteration via row augmentation is no problem (the row augmented part would also stay the same in each iteration and the Ws would be taken to be one for the augmented part) and works out to be exact. As initialization for the variance weights one would take 1/(y+eps) (R uses eps=0.1 in the glm.fit source code in fact) and this also gives estimates that are already be really close to the final estimate in the very first iteration (which is why my weighted nonnegative least squares estimate above works so well - it corresponds to 1 iteration of iteratively reweighted nnls to fit a nonnegative identity link Poisson GLM).

• I think the answer is no; because I don't think it is possible to do ridge regression in the poisson model by row augmentation. The row augmentation in the cited stackexchange question was derived using that the log likelihood is of the form (y-Xb)^T(y-Xb) - something that is not the case for the poisson distribution. I hope someone can give a real proof, though! Apr 28 '19 at 0:45
• Hmm yes I feared that it would not be exact, but the numerical calculations at least show it to be almost identical. So it would still be useful to know if one could say it would be approx. identical and how one should refer to the estimator I suggest. Note that the weighted NNLS estimate above is also almost identical to the nonnegative Poisson estimate (and much better than a regular NNLS estimate) - and that adding a ridge penalty to weighted nnls via row augmentation is mathematically exact. So I could always use that as an approximation to nonnegative Poisson ridge regression as well... Apr 29 '19 at 10:33
• @svendvn Just figured that in the end it may be exact - reason is that you can readily fit nonnegative GLMs using iteratively reweighted nonnegative least squares (just like the regular IRLS algo but substituting the regular solve / lm.wfit by weighted nnls) and for identity link you just have to regress yhat on your covariate matrix and the weights W are just 1/expected variance (ie in the IRLS algo, bwlewis.github.io/GLM, z=yhat and W=1/expected variance). May 15 '19 at 12:05
• @svendvn So given that we are just doing weighted (nonnegative) least squares anyway, adding a ridge penalty in each iteration via row augmentation is no problem (the row augmented part would also stay the same in each iteration and the Ws would be taken to be one for the augmented part) and works out to be exact. May 15 '19 at 12:05
• @svendvn As initialization for the variance weights one would take 1/(y+eps) (R uses eps=0.1) and this also gives estimates that are already be really close to the final estimate in the very first iteration (which is why my weighted nonnegative least squares estimate above works so well - it corresponds to 1 iteration of iteratively reweighted nnls to fit a nonnegative identity link Poisson GLM). May 15 '19 at 12:05