Covariance matrix of estimates using lsqnonneg

Hello,

I want to solve the equation Ax=b for x>=0. How can I derive the covariance matrix of the estimates x. Using ordinary least-squares technique, the covariance matrix is given by

Q = inv(A)
x = Q*b

but how can I get the covariance matrix by using lsqnonneg?

xpos = lsqnonneg(A,b)

/Micha

[Matt J]

I think you mean that the covariance matrix is given by inv(A.'*A).

[John D'Errico]

That is not true for a lsqnonneg problem, IF any of the constraints are active!

[Matt J]

Sure. I meant in the unconstrained case...

Answer 1

[Matt J]

Rewrite the problem with a change of variables

min. f(z)  = || A*z.^2-b||^2

This is now an unconstrained problem, so you can compute the covariance of z in the usual way from the Hessian of f(z) at z=sqrt(abs(xpos)). Then you can map covariance in z to covariance in x.

[John D'Errico]

That gives you only the covariance matrix of z=x^2, not of x. The two will be quite different if any of the constraints were active, thus z was essentially zero! So even if the off-diagonal terms in the resulting covariance matrix of z are sufficiently small that the diagonals of the covariance matrix are valid to be used here, you still have a significant problem.

The point is, suppose that you know the variance of a random variable z, for z very near zero. Thus assume that z has mean mu, with mu essentially zero. You also know the variance of z.

Now compute the variance of sqrt(z). A simple linear approximation will not suffice here, since there is a singularity at zero. In fact, simply assuming the result is normally distributed, with a known variance will be flat out wrong. And if you intend to compute the variance of x, the final goal is always to compute confidence intervals on x. Since the traditional confidence intervals assume a normally distributed random variable x, you have a problem.

So, yes, you can compute a distribution of x, when z is near zero. But it won't be normally distributed.

Effectively, you create this problem by using the change of variables to solve the optimization problem. It still leaves you with no traditionally simple estimate of the uncertainty in x.

You cannot get around the fundamental problem so easily, in that the simple, standard approximations to compute the variance of x are not valid when a constraint is active. Nothing stops you from blithely going ahead and doing that computation. But it may well be a meaningless result.

This is why I suggested some time ago when the question was posed before by this asker, that the best solution is to use a method like the jackknife to compute confidence intervals here. At least if you want valid confidence intervals. If the goal is simply to compute random numbers for confidence intervals, then I suggest either 17 or 42.

[Matt J]

I'm not sure why it would be so hard to extrapolate a meaningful confidence interval from z. I can see that x=sqrt(|z|) will not be normally distributed. But if we're saying that z itself can still be modeled as normally distributed, then we can still deduce from standard confidence intervals that for example | z |<3*sigma with 99.8% probability. From there, we can deduce that 0<x<sqrt(3*sigma) with the same probability. So yes, it's a non-standard confidence interval in x, but does that mean it's not a useful one?

[Matt J]

Michael commented:

Thank you for your suggestions. I found a related question here (lsqnonneg function was not named). The bootstrap method is suggested as an approximation of the uncertainties of the estimates. John reports about reliable results.