# GLM.jl LogisticRegression errors: matrix is not positive definite; Cholesky factorization failed

**URL:** <https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981>\
**Category:** Statistics\
**Tags:** question, glm\
**Created:** [September 20, 2019, 4:37am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981 "2019-09-20T04:37:02Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sijun](https://avatars.discourse-cdn.com/v4/letter/s/b2d939/32.png) [@Sijun](https://discourse.julialang.org/u/Sijun)\
**Post date:** [September 20, 2019, 4:37am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/1 "2019-09-20T04:37:02Z")

</div>

The following errors:

```julia
using DataFrames, GLM

df = DataFrame(x1=[1,2,3,4], x2=[1,2,3,4], y=[1,1,0,0])

mdl = glm(@formula(y~x1+x2), df, Binomial(), LogitLink())

predict(mdl, df[:, [:x1, :x2]])

```

> ERROR: LoadError: PosDefException: matrix is not positive definite; Cholesky factorization failed.

On the other hand, logistic regression is successfully done with ScikitLearn:

```julia

using DataFrames, ScikitLearn

@sk_import linear_model: LogisticRegression

df = DataFrame(x1=[1,2,3,4], x2=[1,2,3,4], y=[1,1,0,0])

mdl = fit!(LogisticRegression(), Matrix(df[:, [:x1, :x2]]), df.y)

ps = predict_proba(mdl, Matrix(df[:, [:x1, :x2]]))

```

Multi-collinearity cannnot be an issue for logistic regression unlike linear regression.

Then, any idea why GLM fails in this case?

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [September 20, 2019, 5:16am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/2 "2019-09-20T05:16:43Z")

</div>

Yes, collinearity can be a problem for nonlinear models, including logistic regression. The coefficients are simply not identified in your case, because x2 is an exact copy of x1. If ScikitLearn gives an answer, it must be eliminating either x1 or x2, or doing something similar to deal with the problem. If the index function is t(x1_a+x2_b), you are probably getting an estimate of a+b, which may be misleading if you don’t notice the collinearity problem. Or, it might restrict a=b, which will also give a biased estimate. In cases like this, I prefer what GLM does, an error message that helps you find the problem.

---

<div class="post-metadata">

**Author:** ![xiaodai](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xiaodai/32/15937_2.png) [@xiaodai](https://discourse.julialang.org/u/xiaodai)\
**Post date:** [September 20, 2019, 6:29am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/3 "2019-09-20T06:29:43Z")

</div>

Scikit learn glm isn’t vanilla glm. It is glm with L2 penalty with magic number 1 as the hyper parameter

---

<div class="post-metadata">

**Author:** ![Sijun](https://avatars.discourse-cdn.com/v4/letter/s/b2d939/32.png) [@Sijun](https://discourse.julialang.org/u/Sijun)\
**Post date:** [September 20, 2019, 7:02am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/4 "2019-09-20T07:02:54Z")

</div>

Linear regression entails matrix inversion, and this is the mechanism via which collinearity affects linear regression in case that the matrix is singular.

However in logistic regression, the estimation of coefficients is based on some likelihood function instead of normal equation as is used in linear regression.

Since there is no closed form expression for obtaining the solution, iterative calculation is instead used for the estimation. There is no way collinearity kicks in the process.

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [September 20, 2019, 7:16am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/5 "2019-09-20T07:16:31Z")

</div>

I think you need [https://github.com/JuliaStats/GLM.jl/pull/314](https://github.com/JuliaStats/GLM.jl/pull/314).

---

<div class="post-metadata">

**Author:** ![Sijun](https://avatars.discourse-cdn.com/v4/letter/s/b2d939/32.png) [@Sijun](https://discourse.julialang.org/u/Sijun)\
**Post date:** [September 20, 2019, 7:22am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/6 "2019-09-20T07:22:36Z")

</div>

Thanks for the information. But it is still unclear how logistic regression can suffer from collinearity unless somehow GLM’s logistic regression algorithm employs matrix calculation like in linear regression( It seems because Cholesky decomposition was mentioned in the error message)

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [September 20, 2019, 7:30am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/7 "2019-09-20T07:30:15Z")

</div>

This is mistaken. Estimation of nonlinear models does entail solving first order conditions, numerically, rather than analytically, as can be done for linear models. With perfect collinearity, there is more than one solution to the first order conditions, and the matrix of second derivatives is singular. Most iterative gradient-based optimization methods (e.g., BFGS) will not work. It may be the case that certain procedures can still produce an estimate (e.g., regularization methods, as noted in another response, or using a derivative free optimization method) but the estimate will not be consistent, in general.

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [September 20, 2019, 7:33am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/8 "2019-09-20T07:33:01Z")

</div>

The algorithm for fitting generalized linear models in the `GLM` package is iteratively reweighted least squares (IRLS) which performs a weighted linear least squares fit at each iteration. Thus a rank-deficient model matrix means that the coefficients are unidentifiable in a GLM as well as in an LM.

The identifiability is a property of the linear predictor, Xβ, where X is the model matrix and β is the coefficient vector. It doesn’t really matter whether the linear predictor is also the fitted value, as in a linear model, or if an inverse link function is applied to it, as in logistic regression.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [September 20, 2019, 9:35am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/9 "2019-09-20T09:35:52Z")

</div>

Possibly useful: [GitHub - timholy/PositiveFactorizations.jl: Positive-definite "approximations" to matrices](https://github.com/timholy/PositiveFactorizations.jl)

---

<div class="post-metadata">

**Author:** ![spinkney](https://avatars.discourse-cdn.com/v4/letter/s/dec6dc/32.png) [@spinkney](https://discourse.julialang.org/u/spinkney)\
**Post date:** [September 20, 2019, 2:31pm UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/10 "2019-09-20T14:31:08Z")

</div>

This is a great answer. For such situations - with the consistency caveat applied - I’ve applied the [Moore-Penrose Inverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse) with success…

---

<div class="post-metadata">

**Author:** ![Sijun](https://avatars.discourse-cdn.com/v4/letter/s/b2d939/32.png) [@Sijun](https://discourse.julialang.org/u/Sijun)\
**Post date:** [September 20, 2019, 2:59pm UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/11 "2019-09-20T14:59:21Z")

</div>

Thanks a lot. It would be useful if GLM could produce p-value and confidence interval along with model coefficients rather than aborts in case there is a collinearity (like in R) so that one can drop collinear predictors.

---

<div class="post-metadata">

**Author:** ![Shuya\_Fan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuya_fan/32/37008_2.png) [@Shuya\_Fan](https://discourse.julialang.org/u/Shuya_Fan)\
**Post date:** [June 8, 2022, 9:10am UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/12 "2022-06-08T09:10:35Z")

</div>

I have the same error when I ran the GLMER model: `dfAllform = @formula(ys ~ 1 + bio_pd+ (1|species)+ (1|region)) gm1 = fit(MixedModel, dfAllform, dfAll, Binomial(), CloglogLink(), fast = true)`  
In my model, there is only one fixed effect, so collinearity won’t be the cause of the error. Does anyone know there are other possible reasons?

---

<div class="post-metadata">

**Author:** ![dmbates](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmbates/32/44_2.png) [@dmbates](https://discourse.julialang.org/u/dmbates)\
**Post date:** [June 8, 2022, 2:12pm UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/13 "2022-06-08T14:12:28Z")

</div>

A similar problem reported in [LoadError: PosDefException: matrix is not positive definite; Cholesky factorization failed · Issue #611 · JuliaStats/MixedModels.jl · GitHub](https://github.com/JuliaStats/MixedModels.jl/issues/611) turned out to be due to scaling of the response and model matrix - because of weights in that case. That example was a linear mixed model (LMM). With generalized linear mixed models (GLMMs) the problem is even more likely because of the iteratively re-weighted least squares (IRLS) algorithm and its penalized version (PIRLS). Obtaining reasonable starting estimates to avoid this is not easy.

---

<div class="post-metadata">

**Author:** ![Shuya\_Fan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shuya_fan/32/37008_2.png) [@Shuya\_Fan](https://discourse.julialang.org/u/Shuya_Fan)\
**Post date:** [June 9, 2022, 3:46pm UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/14 "2022-06-09T15:46:58Z")

</div>

Thanks for your kindly reply. Actually, I still don’t quite understand what should I do to solve this problem. Accroding to your suggestion, the response variable (`ys` in my model) should not be scaled or the parameters should be set to certain initial values, then I can deal with this error by the two methods? I am sorry that I barely know the specific fitting process of the GLMM model. If I misunderstood, please let me know. Thanks again.

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [June 9, 2022, 10:04pm UTC](https://discourse.julialang.org/t/glm-jl-logisticregression-errors-matrix-is-not-positive-definite-cholesky-factorization-failed/28981/15 "2022-06-09T22:04:16Z")

</div>

It’s best if you open a new thread because this is getting further and further away from the original GLM.jl issue. The problem with fitting GLMMs is broadly similar to the problem in fitting GLMs, but there are a lot of differences in the details.

It’s also very hard for us to further debug your problem without having access to the data to try to reproduce the problem. If you are unable to publicly share your data, but you could share the data with me or @dmbates, you can contact us via direct message and provide a private link to the data.