GhostBasil very slow for p=1140 problem

[biona001]

Hi James,

Here is an example where GhostBasil runs slowly. This region has 190 variables, and I generate 5 knockoff copies, so overall A is 1140*1140. On sherlock with 12 cores, it takes ~200 seconds to converge.

Could you check if the output below matches your expectation? Hopefully I didn't misuse your package in some way.

• Ci.txt
• lambdas.txt
• r.txt
• S.txt

library(Matrix)
library(ghostbasil)

# read data
Ci <- as.matrix(read.table('Ci.txt'))
S <- as.matrix(read.table('S.txt'))
r <- scan('r.txt')
lambdas <- scan('lambdas.txt')

# form ghostbasil inputs
S <- BlockMatrix(list(S))            # dim(S) = 190 by 190
A <- BlockGroupGhostMatrix(Ci, S, 6) # dim(A) = 1140 by 1140

# run ghostbasil
result <- ghostbasil(A, r, delta.strong.size=500, user.lambdas=lambdas,
    max.strong.size = nrow(A), n.threads=12, use.strong.rule=F)

Some observations:

• The model seems to converged (no error message printed and result$error is "")
• result$rsqs returns a vector of NaN

  result$rsqs
  [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  [19] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  [37] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  [55] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  [73] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  [91] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

• The number of non-zero beta is not uniformly increasing

  colSums(result$betas != 0)
  [1]   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  [19]   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  [37]   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  [55]   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
  [73]   0   0   0   0   0   0   0   0   0   0   0   0   0   1   8  11  18  31
  [91]  37  48  70 104 147 195 286 428 706   5

[biona001]

A few more observations

• Since my block size is 1, there is nothing to parallelize here, so I should actually just use n.threads=1
• I tried delta.strong.size=10 and the same problem is still running after one hour

[JamesYang007]

Damn I'm inclined to believe there's a 🐛. Thanks for finding this example!! I'll take a look asap. I probably won't be able to get to it till next week though.

[JamesYang007]

Ok started looking into this! I don't think it's a bug anymore. Here are some of my findings:

• S and Ci are not symmetric (close to symmetric though). I don't know if this matters so much, but for numerical stability, it might be worth symmetrizing: S = (S + t(S)) / 2 (similarly for Ci).
• MORE IMPORTANTLY: the resulting A is not PSD!

  A.dense <- A$to_dense()
  A.dense.eig = eigen(A.dense, symmetric=T)
  tail(A.dense.eig$values, 10)

  [1] -0.01602446 -0.02566999 -0.02675701 -0.02717474 -0.02850788 -0.02912589
  [7] -0.03217679 -0.03340235 -0.03437624 -0.03684428

  I get very similar results if I symmetrize A.dense directly. These eigenvalues are significantly negative. I did not change your code snippet - just added these three lines before calling ghostbasil.

• I noticed that your lambda sequence is increasing. It should be decreasing (as mentioned in the email thread as well)!
• I modified the A such that I truncate the negative eigenvalues to 0 and then I add 1e-10 to the diagonal to make it slightly PD:

  dd = sapply(A.dense.eig$values, function(x) max(x, 0)) + 1e-10
  A.dense.trunc = A.dense.eig$vectors %*% (dd * t(A.dense.eig$vectors))

  Now, running ghostbasil gives good result:

  result.dense <- ghostbasil(A.dense.trunc, r, alpha=1, delta.strong.size=500, user.lambdas=lambdas,
  max.strong.size = nrow(A.dense.trunc), n.threads=1, use.strong.rule=T)

  For example, here's the (relative) rsq:

  [1] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [6] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [11] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [16] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [21] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [26] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [31] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [36] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [41] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [46] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [51] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [56] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [61] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [66] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [71] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [76] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [81] 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
  [86] 0.0003982219 0.0044517538 0.0110897322 0.0204761087 0.0345857061
  [91] 0.0526677093 0.0725795312 0.0970495241 0.1328570606 0.1861651392
  [96] 0.2550997130 0.3586244934 0.5188728742 0.7549665432 1.0000000000

I suspect that something's weird in the construction of Ci or S because the group-knockoff construction should ensure that the knockoff matrix A is PSD at the very least. When you solved for S, were you assuming only 1 knockoff or something?

[biona001]

Hi James, thanks for taking a look. A definitely should be PSD, so I must have messed up somewhere, let's see...

[biona001]

After reversing lambda and making sure A is PSD, the same problem now runs <1 second. Thanks for the insights, and sorry for the troubles.