Error in solve.default(comp$XX): system is computationally singular

[kfuku52]

Sorry to keep bothering you. I found another error which happens to some of my relatively large datasets (typically >150 leaves). Like the problem in #29, the error can be circumvented by specifying fixed.alpha=TRUE in estimate_convergent_regimes(), although the error message is different. I would appreciate any comments or suggestion.

Input data: traits.txt tree.txt

Error message:

Error in solve.default(comp$XX): system is computationally singular: reciprocal condition number = 4.72733e-22
Traceback:

1. estimate_convergent_regimes(fit_ind, criterion = "AICc", method = "backward", 
 .     fixed.alpha = FALSE, nCores = 1)
2. estimate_convergent_regimes_surface(model, opt = opt)
3. cmp_model_score_CR(tr, Y, regimes, model$alpha, opt = opt)
4. cmp_AICc_CR(tree, Y, conv.regimes = regimes, alpha = alpha, opt = opt)
5. phylolm_interface_CR(tree, matrix(Y[, i]), conv.regimes, alpha = alpha[[i]], 
 .     opt = opt)
6. phylolm_CR(Y ~ preds - 1, phy = tr, model = "OUfixedRoot", sc = shift.configuration, 
 .     cr = conv.regimes, starting.value = alpha, upper.bound = opt$alpha.upper.bound, 
 .     lower.bound = alpha/100)
7. optim(logstart, fn = minus2llh, method = "L-BFGS-B", lower = loglower, 
 .     upper = logupper, ...)
8. (function (par) 
 . fn(par, ...))(-22.855602561775)
9. fn(par, ...)
10. loglik(prm, y, X = preds)
11. solve(comp$XX)
12. solve.default(comp$XX)

Code to reproduce:

library(l1ou)
tree = read.tree("tree.txt")
trait_matrix = read.table("traits.txt", sep="\t")
adj_data = adjust_data(tree=tree, Y=trait_matrix, normalize = FALSE)
fit_ind = estimate_shift_configuration(tree=adj_data$tree, Y=adj_data$Y, criterion="AICc", root.model="OUrandomRoot", nCores=4, rescale=TRUE)

# no error
fit_conv = estimate_convergent_regimes(fit_ind, criterion="AICc", method="backward", fixed.alpha=TRUE, nCores=4)

# error
fit_conv = estimate_convergent_regimes(fit_ind, criterion="AICc", method="backward", fixed.alpha=FALSE, nCores=4)

Many thanks.

Kenji

[khabbazian]

Hi Kenji,

Thank you for sending the error. I know what causes the error. I need some time to figure out how to fix it. So

If it is okay I will work on it more over the weekend and get back to you afterward.

Best regards,

Mohammad

---

From: Kenji Fukushima notifications@github.com Sent: Wednesday, May 31, 2017 1:23:43 PM To: khabbazian/l1ou Cc: Subscribed Subject: [khabbazian/l1ou] Error in solve.default(comp$XX): system is computationally singular (#30)

Sorry to keep bothering you. I found another error which happens to some of my relatively large datasets (typically >150 leaves). Like the problem in #29https://github.com/khabbazian/l1ou/issues/29, the error can be circumvented by specifying fixed.alpha=TRUE in estimate_convergent_regimes(), although the error message is different. I would appreciate any comments or suggestion.

Input data: traits.txthttps://github.com/khabbazian/l1ou/files/1039445/traits.txt tree.txthttps://github.com/khabbazian/l1ou/files/1039446/tree.txt

Error message:

Error in solve.default(comp$XX): system is computationally singular: reciprocal condition number = 4.72733e-22 Traceback:

1. estimate_convergent_regimes(fit_ind, criterion = "AICc", method = "backward", . fixed.alpha = FALSE, nCores = 1)
2. estimate_convergent_regimes_surface(model, opt = opt)
3. cmp_model_score_CR(tr, Y, regimes, model$alpha, opt = opt)
4. cmp_AICc_CR(tree, Y, conv.regimes = regimes, alpha = alpha, opt = opt)
5. phylolm_interface_CR(tree, matrix(Y[, i]), conv.regimes, alpha = alpha[[i]], . opt = opt)
6. phylolm_CR(Y ~ preds - 1, phy = tr, model = "OUfixedRoot", sc = shift.configuration, . cr = conv.regimes, starting.value = alpha, upper.bound = opt$alpha.upper.bound, . lower.bound = alpha/100)
7. optim(logstart, fn = minus2llh, method = "L-BFGS-B", lower = loglower, . upper = logupper, ...)
8. (function (par) . fn(par, ...))(-22.855602561775)
9. fn(par, ...)
10. loglik(prm, y, X = preds)
11. solve(comp$XX)
12. solve.default(comp$XX)

Code to reproduce:

library(l1ou) tree = read.tree("tree.txt") trait_matrix = read.table("traits.txt", sep="\t") adj_data = adjust_data(tree=tree, Y=trait_matrix, normalize = FALSE) fit_ind = estimate_shift_configuration(tree=adj_data$tree, Y=adj_data$Y, criterion="AICc", root.model="OUrandomRoot", nCores=4, rescale=TRUE)

no error

fit_conv = estimate_convergent_regimes(fit_ind, criterion="AICc", method="backward", fixed.alpha=TRUE, nCores=4)

error

fit_conv = estimate_convergent_regimes(fit_ind, criterion="AICc", method="backward", fixed.alpha=FALSE, nCores=4)

Many thanks.

Kenji

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/khabbazian/l1ou/issues/30, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AH3biz2QiwLUuKG1DZzwKglu9qcy62yHks5r_aIfgaJpZM4Nr8U_.

[kfuku52]

Thank you! I look forward to hearing your feedback. -Kenji

[khabbazian]

I updated the package with the fix. I couldn't run your code to make sure. Could you let me know if you still get the same error?

---

From: Kenji Fukushima notifications@github.com Sent: Thursday, June 1, 2017 10:50:20 AM To: khabbazian/l1ou Cc: Mohammad KHABBAZIAN; Comment Subject: Re: [khabbazian/l1ou] Error in solve.default(comp$XX): system is computationally singular (#30)

Thank you! I look forward to hearing your feedback. -Kenji

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/khabbazian/l1ou/issues/30#issuecomment-305517332, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AH3bi0rnIuyfKmtqbK58YWT34H1TMYIHks5r_s-sgaJpZM4Nr8U_.

[kfuku52]

Thank you. I updated the package and tried again but got another error in estimate_convergent_regimes().

Error in optim(logstart, fn = minus2llh, method = "L-BFGS-B", lower = loglower,  : 
  L-BFGS-B needs finite values of 'fn'
In addition: Warning message:
In mclapply(run.list, FUN = function(X) { :
  all scheduled cores encountered errors in user code

[khabbazian]

Could you please update the package and run your code again?

The problem is while optim tries to optimize the loglik, for the small values of alpha the design matrix becomes singular ( the condition number get less thatn .Machine$double.eps). I changed the tol in the solve(.). Let see if that helps.

---

From: Kenji Fukushima notifications@github.com Sent: Tuesday, June 6, 2017 6:19:55 PM To: khabbazian/l1ou Cc: Mohammad KHABBAZIAN; Comment Subject: Re: [khabbazian/l1ou] Error in solve.default(comp$XX): system is computationally singular (#30)

Thank you. I updated the package and tried again but got another error in estimate_convergent_regimes().

Error in optim(logstart, fn = minus2llh, method = "L-BFGS-B", lower = loglower, : L-BFGS-B needs finite values of 'fn' In addition: Warning message: In mclapply(run.list, FUN = function(X) { : all scheduled cores encountered errors in user code

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/khabbazian/l1ou/issues/30#issuecomment-306633903, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AH3bixsNg0UAPDAlXNojw784xiDFirlmks5sBdCLgaJpZM4Nr8U_.

[cecileane]

Another idea, Mohammad, would be to call the log-likelihood function within a try or tryCatch block, and return a log-likelihood of -Inf if the design matrix becomes singular. That would force the estimate of α to be reasonably not too small.

[khabbazian]

I did something like that before this one but we got the following error

Error in optim(logstart, fn = minus2llh, method = "L-BFGS-B", lower = loglower, :

L-BFGS-B needs finite values of 'fn'

I will check if we get similar error if it returns NA.

---

From: Cécile Ané notifications@github.com Sent: Friday, June 9, 2017 9:46:20 AM To: khabbazian/l1ou Cc: Mohammad KHABBAZIAN; Comment Subject: Re: [khabbazian/l1ou] Error in solve.default(comp$XX): system is computationally singular (#30)

Another idea, Mohammad, would be to call the log-likelihood function within a try or tryCatch block, and return a log-likelihood of -Inf if the design matrix becomes singular. That would force the estimate of α to be reasonably not too small.

― You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/khabbazian/l1ou/issues/30#issuecomment-307392950, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AH3bi181h0IG2JDM8aX5o5e8KHTmJ8JRks5sCUysgaJpZM4Nr8U_.

[kfuku52]

The updated version of estimate_convergent_regimes() is still running (>6 hrs from start, 16 cores allocated). I will get back to you when it finishes.

[khabbazian]

I guess because I changed the tol parameters of solve from 1e-16 to 1e-64 It became slower. Look forward.

---

From: Kenji Fukushima notifications@github.com Sent: Friday, June 9, 2017 6:41:38 PM To: khabbazian/l1ou Cc: Mohammad KHABBAZIAN; Comment Subject: Re: [khabbazian/l1ou] Error in solve.default(comp$XX): system is computationally singular (#30)

The updated version of estimate_convergent_regimes() is still running (>6 hrs from start, 16 cores allocated). I will get back to you when it finishes.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/khabbazian/l1ou/issues/30#issuecomment-307516258, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AH3biyZMDTE6rIZWZRl5VJb-FVbht9sKks5sCcoigaJpZM4Nr8U_.

[kfuku52]

The process finished successfully without error. Thank you very much! Also, I would appreciate it if you could help fix #29 when you have time.

[khabbazian]

Thank you. I will work on that too.

---

From: Kenji Fukushima notifications@github.com Sent: Sunday, June 11, 2017 11:49:12 AM To: khabbazian/l1ou Cc: Mohammad KHABBAZIAN; Comment Subject: Re: [khabbazian/l1ou] Error in solve.default(comp$XX): system is computationally singular (#30)

The process finished successfully without error. Thank you very much! Also, I would appreciate it if you could help fix #29https://github.com/khabbazian/l1ou/issues/29 when you have time.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/khabbazian/l1ou/issues/30#issuecomment-307637938, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AH3biwJ499AwvK6IgwFO5eCnIrj5Eg8Qks5sDAx4gaJpZM4Nr8U_.