GSL ERROR: matrix is singular in lu.c at line 266 errno 1

[pjotrp]

On Fri, Sep 21, 2018 at 03:58:45AM -0700, Xiang Zhou wrote:

Have you tried the early 0.94.1 version ([1]http://xzlab.org/software/
gemma-0.94.1/gemma)? For multiple trait analysis, I encountered a few
settings before where the early version seems to work slightly more
stable.

As per your suggestion I ran 0.94, 0.96 and 0.98 with the correlated phenotypes using

~/opt/gemma-0.94/gemma-v0.94 -p data/correlated_phenotypes/Ysim_reg_gemma.txt     -g data/correlated_phenotypes/Genotypes_gemma.csv -d data/c
orrelated_phenotypes/Kinship_eigenval_gemma.txt     -u data/correlated_phenotypes/Kinship_eigenvec_gemma.txt -lmm 2 -n 1 9 4 6 10 -o corrpheno9
4

and there is no real difference between 0.94 and 0.96. They both give a result.

0.98, however, gives with the -no-check option

  ERROR: Enforce failed for LU determinant is zero -> LU is not invertable in src/lapack.cpp at line 339 in LULndet

Now we could disable that check too, but I think it is correct and
should stop GEMMA. Removing that check we get the familiar

  GSL ERROR: matrix is singular in lu.c at line 266 errno 1

which points out there is a problem LUInvert in the mvlmm::CalcDev function

  // Invert Hessian.
  int sig;
  gsl_permutation *pmt = gsl_permutation_alloc(v_size * 2);

  LUDecomp(Hessian, pmt, &sig);
  LUInvert(Hessian, pmt, Hessian_inv);

This is the same code you can find in v0.96 at

https://github.com/genetics-statistics/GEMMA/blob/v0.96/src/mvlmm.cpp#L2401

Compiling with GSLv1 gives the same error as GSLv2.

Why the older versions pass this check is (at this point) a mystery to me.

Hannah, can you try and see if there is a difference between 0.94 and
0.96 works for all your cases? It would be great if we find one that
works in one of those versions and fails in the other.

I am close to putting out a GEMMA 0.98 release, but I think we ought
to fix this issue.

[pjotrp]

@HannahVMeyer an you try and see if there is a difference between 0.94 and
0.96 works for all your cases? It would be great if we find one that
works in one of those versions and fails in the other.

[pjotrp]

OK, I found something of interest. Gemma 0.96 when built on my system fails, but the binary 0.96 release provided by @pcarbo passes also for all 10 phenotypes. Peter wrote

Executable gemma.linux was built with an Intel Xeon E5-2680v4 ("Broadwell") 2.4GHz processor, Scientific Linux 7 (64-bit) and gcc 4.8.5, and statically linked to glibc 2.14, atlas 3.10.3 and gsl 1.16.

So, the difference here is atlas and a much older GSL. The actual libs imported are

ldd ~/opt/gemma-0.96/gemma.linux 
        linux-vdso.so.1 (0x00007ffdd57ba000)
        libz.so.1 => /lib/x86_64-linux-gnu/libz.so.1 (0x00007f11b0d0f000)
        libgfortran.so.3 => /usr/lib/x86_64-linux-gnu/libgfortran.so.3 (0x00007f11b09f1000)
        libstdc++.so.6 => /usr/lib/x86_64-linux-gnu/libstdc++.so.6 (0x00007f11b066f000)
        libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007f11b036b000)
        libgcc_s.so.1 => /lib/x86_64-linux-gnu/libgcc_s.so.1 (0x00007f11b0154000)
        libpthread.so.0 => /lib/x86_64-linux-gnu/libpthread.so.0 (0x00007f11aff37000)
        libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f11afb99000)
        libquadmath.so.0 => /usr/lib/x86_64-linux-gnu/libquadmath.so.0 (0x00007f11af95a000)
        /lib64/ld-linux-x86-64.so.2 (0x00007f11b0f29000)

Also of interest is that 0.94 is faster than 0.96. At least for this dataset.

[pjotrp]

Good news, I replicated the problem wich compares with v94 output. v96 with GSL1.16 returns:

env GUIX_PACKAGE_PATH=~/izip/git/opensource/genenetwork/guix-bioinformatics/ ~/.config/guix/current/bin/guix environment -C guix --ad-hoc gcc gdb gfortran:lib gsl1 eigen atlas lapack zlib bash ld-wrapper perl --substitute-urls="https://berlin.guixsd.org  https://mirror.hydra.gnu.org http://guix.genenetwork.org"

time ./bin/gemma -p ~/tmp/hannah/Ysim_reg_gemma.txt     -g ~/tmp/hannah/Genotypes_gemma.csv -d ~/tmp/hannah/Kinship_eigenval_gemma.txt     -u ~/tmp/hannah/Kinship_eigenvec_gemma.txt -lmm 2 -n 1 9 4 6 10 -o corrpheno96
Reading Files ... 
## number of total individuals = 1000
## number of analyzed individuals = 1000
## number of covariates = 1
## number of phenotypes = 5
## number of total SNPs = 100
## number of analyzed SNPs = 100
REMLE estimate for Vg in the null model: 
0.8091
0.1435  0.0660
-0.1456 -0.0664 0.1165
0.0907  0.0362  -0.0830 0.0639
-0.5533 -0.0985 0.1171  -0.0765 0.4131
se(Vg): 
0.0736
0.0355  0.0299
0.0421  0.0264  0.0422
0.0444  0.0279  0.0430  0.0479
0.0516  0.0254  0.0295  0.0313  0.0389
REMLE estimate for Ve in the null model: 
0.0334
0.0267  0.4502
-0.0529 -0.0883 0.6004
0.0683  0.1413  -0.6502 0.8400
-0.0128 -0.0041 -0.0706 0.0880  0.0419
se(Ve): 
0.0378
0.0252  0.0329
0.0299  0.0281  0.0456
0.0332  0.0319  0.0494  0.0583
0.0266  0.0186  0.0215  0.0242  0.0208
REMLE likelihood = -3813.3258
MLE estimate for Vg in the null model: 
0.8082
0.1434  0.0660
-0.1455 -0.0664 0.1165
0.0906  0.0362  -0.0830 0.0639
-0.5528 -0.0984 0.1170  -0.0764 0.4127
se(Vg): 
0.0735
0.0355  0.0298
0.0420  0.0264  0.0421
0.0443  0.0278  0.0430  0.0478
0.0515  0.0254  0.0295  0.0313  0.0388
MLE estimate for Ve in the null model: 
0.0334
0.0267  0.4497
-0.0529 -0.0882 0.5997
0.0682  0.1412  -0.6495 0.8391
-0.0128 -0.0040 -0.0705 0.0879  0.0418
se(Ve): 
0.0378
0.0252  0.0329
0.0298  0.0281  0.0455
0.0331  0.0318  0.0493  0.0582
0.0265  0.0185  0.0215  0.0242  0.0207
MLE likelihood = -3814.3913
Reading SNPs  ==================================================100.00%

real    0m1.486s
user    0m1.448s
sys     0m0.036s

So, again, it is a GSL thing.

[pjotrp]

v96 with openblas it gives the same result. Great, now I can do a git bisect.

[pjotrp]

9952786 is the culprit.

[ohines]

Hi,
I appear to have encountered a the same problem using v98.1
below is my function call and read out
thank you in advance

./gemma -g input/1000G_P3V5_PID_chr21.dose.bimbam -p input/pheno_final.txt -c input/cov_final.txt -k output/GSM_K_Final.sXX.txt -lmm 4 -n 1 2 3 -o results
GEMMA 0.98.1 (2018-12-10) by Xiang Zhou and team (C) 2012-2018
Reading Files ...
## number of total individuals = 5654
## number of analyzed individuals = 3020
## number of covariates = 6
## number of phenotypes = 3
## number of total SNPs/var        =   653791
## number of analyzed SNPs         =   137887
Start Eigen-Decomposition...
GSL ERROR: matrix is singular in lu.c at line 266 errno 1

[pjotrp]

Did you search the issue tracker? Try with less covariates. Usually it is because they are correlated.

[pjotrp]

https://github.com/genetics-statistics/GEMMA/issues?utf8=%E2%9C%93&q=is%3Aissue+singular

[ohines]

Hi Pjotr,
Thanks for getting back to me so quickly.
I reran the command with no covariate file, but recovered the same issue.

Reading Files ...
## number of total individuals = 5654
## number of analyzed individuals = 3020
## number of covariates = 1
## number of phenotypes = 3
## number of total SNPs/var        =   653791
## number of analyzed SNPs         =   138110
Start Eigen-Decomposition...
GSL ERROR: matrix is singular in lu.c at line 266 errno 1

Then again with no covariate file and only phenotype1, which was successful

GEMMA 0.98.1 (2018-12-10) by Xiang Zhou and team (C) 2012-2018
Reading Files ...
## number of total individuals = 5654
## number of analyzed individuals = 3020
## number of covariates = 1
## number of phenotypes = 1
## number of total SNPs/var        =   653791
## number of analyzed SNPs         =   138110
Start Eigen-Decomposition...
pve estimate =0.623785
se(pve) =0.030895
================================================== 100%
**** INFO: Done.

The command was also successful with the covariate file included and only 1 phenotype.
Could the issue be that the phenotypes are linearly dependent?
In this case the phenotypes (intervals in an ECG trace) are constructed such that
Phenotype1 = Phenotype2 + Phenotype3
I had thought that each phenotype would be analysed separately

[pjotrp]

Avoid correlated phenotypes/covariates.

[asafpr]

Sorry to dig up an old issue. I ran into the same problem so I used svd and used the U matrix as input so all the phenotypes are uncorrelated and I still got the error, it worked though when I used a subset of the matrix, even 8 out of 9 phenotypes worked. Any idea why? Thanks

[pjotrp]

Not really, sorry.