Numerical instability with general X, D

[huisaddison]

There are numerical instability issues in the general X, D setting which can lead genlasso to terminate unexpectedly early.

See the MRE here: https://gist.github.com/huisaddison/31ea593ae945bb2494033abcf0ddbf1a

Running this MRE yields the output, when svd=TRUE:

1. lambda=107324.521, adding coordinate 146, |B|=1...
2. lambda=56874.590, adding coordinate 1, |B|=2...
3. lambda=0.771, adding coordinate 145, |B|=3...
4. lambda=0.029, adding coordinate 68, |B|=4...
5. lambda=0.024, adding coordinate 69, |B|=5...
6. lambda=0.024, deleting coordinate 68, |B|=4...
7. lambda=0.022, adding coordinate 70, |B|=5...
8. lambda=0.022, deleting coordinate 69, |B|=4...
9. lambda=0.020, adding coordinate 71, |B|=5...
10. lambda=0.019, deleting coordinate 70, |B|=4...
Reached the maximum number of steps (10), skipping the rest of the path.

but when svd=FALSE (the default):

1. lambda=20415224.178, adding coordinate 146, |B|=1...
2. lambda=2604.649, adding coordinate 145, |B|=2...
3. lambda=0.000, deleting coordinate 146, |B|=1...

Sorry for not providing much more information aside from this observation, I haven't dug deep into this myself to fully understand why this is occurring. Some notes:

• Numerical instability introduced when transforming the generalized lasso problem into the equivalent signal approximator problem by multiplying D by X-pinv (despite adding a sizable diagonal ridge). See code:

  genlasso/R/genlasso.R

  Lines 69 to 72 in 5adbeb9

  	x = svd(rbind(X,diag(sqrt(eps),p)))
  	y2 = as.numeric(x$u %*% t(x$u) %*% c(y,rep(0,p)))
  	Xi = x$v %*% (t(x$u) / x$d)
  	D2 = D %*% Xi

• FWIW, in another example (not included, but I can produce one if it's helpful) -- the problem terminates early when svd=TRUE, but runs for longer when svd=FALSE. (And the only difference between these problems is numerical error between the D matrices). See the R Markdown notebook here: https://gist.github.com/huisaddison/da568096fbf01f780e7f06aa9616b8de So the svd=TRUE routine isn't necessarily a "gold standard" for bypassing the numerical issues -- there's something that's going on that I don't fully understand.

[ryantibs]

Thank you @huisaddison for pointing this out ... scary and subtle!

My first thought is that initial linear system defined by D D^T u = D y, where D is the modified analysis operator (defined by the pseudoinverse of your X = P and the original D) must be terribly conditioned. Otherwise I don’t understand how the paths could be so divergent. I think basically different ways of solving this system interact (SVD versus QR) with different versions of the forming “original D” to produce very different solutions.

So what I would do to debug this genlasso issue would just be to get all those modified D’s in memory, compute their condition numbers, and and then try computing all the solutions (SVD versus QR). And then look at them (say, plot them) to see how divergent they are.

(Not suggesting that you should do this ... if @statsmaths or I have free cycles we could do this too ... just recording my own thoughts so I don't forget them ...)

[huisaddison]

A minimal set of code for constructing the modified D's: https://gist.github.com/huisaddison/fe089e9dce7d61dcaeca84e6553664fd

The condition numbers are pretty large, but the matrices look relatively close to each other, especially compared to their absolute magnitudes:

> rcond(modD_ryan)
[1] 4.773753e-19
> rcond(modD_dspline)
[1] 7.59923e-19
> rcond(modD_genlasso)
[1] 4.305832e-19
> max(abs(modD_ryan - modD_dspline))
[1] 1.525879e-05
> max(abs(modD_ryan - modD_genlasso))
[1] 1.525879e-05
> max(abs(modD_dspline - modD_genlasso))
[1] 7.629395e-06
> summary(as.vector(modD_dspline))
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-6.082e+10  0.000e+00  0.000e+00  0.000e+00  0.000e+00  6.289e+10
> summary(as.vector(modD_ryan))
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-6.082e+10  0.000e+00  0.000e+00  0.000e+00  0.000e+00  6.289e+10
> summary(as.vector(modD_dspline))
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-6.082e+10  0.000e+00  0.000e+00  0.000e+00  0.000e+00  6.289e+10
> summary(as.vector(modD_genlasso))
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max.
-6.082e+10  0.000e+00  0.000e+00  0.000e+00  0.000e+00  6.289e+10

I also looked at the modified y's -- negligible differences there.

> max(abs(y2_dspline - y2_genlasso))
[1] 0
> max(abs(y2_ryan - y2_genlasso))
[1] 0
> max(abs(y2_ryan - y2_dspline))
[1] 0

Next steps would be to look at the solution paths themselves using dualpath and dualpathSvd