ARPACK randomly not converging during hwe_normalized_pca

I’m getting

Error summary: IllegalStateException: ARPACK returns non-zero info = 1 Maximum number of iterations taken. (Refer ARPACK user guide for details)

after calling hwe_normalized_pca(mt.GT, k=2) and I have trouble pinpointing the cause of the problem. What’s worse, when I’m calling this function several times in a row on the same input, sometimes it errors out and sometimes it returns without an error. Is that even possible with a proper installation? I have now made myself a jupyter notebook with the following single cell just to prove to myself that I’m not going mad:

mt = hl.read_matrix_table(my_path)
for i in range(12):
            k = 2
        print(f"{i} success")
        print(f"{i} error")

and it prints

0 success
1 success
2 success
3 success
4 error
5 success
6 success
7 success
8 success
9 error
10 success
11 success

where mt is a product of a several filtering steps performed on a genuine vcf.
Also, when I do

    k = 2

on my mt, I always immediately get another kind of error:
Error summary: HailException: Cannot create RowMatrix: filtered entry at row 8111 and col 9,
but I don’t know if it is related to the problem with hwe_normalizd_pca.
Any suggestions would be appreciated:

  1. What qualities of any mt could in general cause the hwe_normalized_pca to not converge?
  2. Should I suspect that there’s something wrong with my installation? I’m running hail locally, on a laptop.
  3. Shape of mt is (18205, 39), its only non-key field is GT, and the col-keys are generic strings, so I’d be happy to share it.

You could make normal PCA work with hl.or_else to set a value for missing entries.

There’s almost certainly no issue with your installation. I can only speculate as to the issues, but I’d guess that your matrix is exhibiting some kind of numerical instability that makes ARPACK’s SVD routines occasionally not converge. Your dataset is quite small. You might try collecting it locally and converting it to a numpy array for use with sklearn.decomposition.PCA. If that succeeds, it’s possible Hail could use whatever parameters sklearn is using to call ARPACK and achieve more consistent results.

Thanks, @danking, for the response!
It’s helpful to hear that reinstalling probably won’t help.
Yup, the dataset is small, partly because I was trying to take a subset of it which still throws the same error as the original :slight_smile:
Not sure where I would put hl.or_else in But, swapping mt.GT for mt.unfilter_entries().GT made the hl.pca return without an error. Similarly

from sklearn.decomposition import PCA
mt_2 = mt.unfilter_entries()
df = mt_2.select_entries(
    n =

returns without an error, but it would complain about NaN without the unfilter_entries (same with say n = hl.or_else(mt_2.GT.n_alt_alleles(), 0)).
Still, unfilter_entries does not help with the hwe_normalized_pca.
Perhaps I’m having a problem with variants that are far from hw equilibrium. E.g. if I filter out low hw p-value before calling hwe_normalized_pca in my jupyter cell:

mt = hl.read_matrix_table("debugging/")
mt = mt.filter_rows(
    hl.agg.hardy_weinberg_test(mt.GT).p_value > 1e-6

it throws away 25 variants and prints

0 success
1 success
2 success
3 error
4 success
5 success
6 success
7 success
8 success
9 success
10 success
11 success

It sounds to me like hardy-weinberg normalization is producing a matrix that is likely to be poorly conditioned causing ARPACK to not converge.

Let’s compare apples to apples:

import hail as hl
mt = hl.read_mtarix_table('debugging/')
mt_for_comparison = hl.methods.pca.hwe_normalize(mt.GT)
# ... run mt_for_comparison through sklearn PCA and hl.pca

How would I import and use the hwe_normalize function? It doesn’t seem to be a part of Hails interface. The following fails

from hail.methods.pca import hwe_normalize
    n = hwe_normalize(mt.GT)

at select_entries with
ExpressionException: Cannot combine expressions from different source objects..

This is not documented because it’s an internal method. This conversation may end up being an impetus for us to document it and expose it. The usage is similar to what Dan did above, but you’d need to import it the way you did. So:

from hail.methods.pca import hwe_normalize

norm_gt = hwe_normalize(mt.GT)

# note it's not a matrix table, its an entry expression. You can export it for sklearn though using: