parallel computing in cv_iht is slow for floating point matrices.

[biona001]

As noted, currently the parallel computing part for the cross validation routine is slower than single core. I don't really know why at the moment, but I suspect it is caused oversubscription with pmap and BLAS mixed together..

For the moment, use parallel = false when you are using single/double-precision matrices.

[kose-y]

I agree that it might be oversubscription with pmap and BLAS. Have you tried running with environment variables JULIA_NUM_THREADS=1, OMP_NUM_THREADS=1, and/or OPENBLAS_NUM_THREADS=1? That would avoid the issue, each process running only a single thread.

[biona001]

The advice I got was to BLAS.set_num_threads(1), but I haven't had time to test it yet. This will get fixed in the next few weeks... hopefully....

[kose-y]

Just in case you are planning to add that line to the main package (under /src), I would recommend not to. That can be an unexpected side effect for some users. I think documenting the issue and using the line in the experiment code in notebooks should be enough.

[biona001]

With recent changes in cross validation, this issue is resolved (it may have been resolved a long time ago but I never checked).

MWE

using MendelIHT, Random, SnpArrays, GLM, LinearAlgebra
BLAS.set_num_threads(1)
Random.seed!(1111)
n = 1000
p = 10000
k = 10
d = Normal
l = IdentityLink()

# simulate genotypes and nongenetic covariate
x = randn(n, p)
z = ones(n, 2) 
z[:, 2] .= rand(0:1, n)
standardize!(@view(z[:, 2:end])) 

# randomly set genetic predictors where causal βᵢ ~ N(0, 1)
true_b = zeros(p) 
true_b[1:k-2] = randn(k-2)
shuffle!(true_b)

# define effect size of non-genetic predictors: intercept & sex
true_c = [1.0; 1.5] 

# simulate phenotype using genetic and nongenetic predictors
prob = GLM.linkinv.(l, x * true_b .+ z * true_c) # note genotype-vector multiplication is done with `xla`
y = [rand(d(i)) for i in prob]

@time mses = cv_iht(y, x, z, path=1:20, verbose=false);
@time mses = cv_iht(y, x, z, path=1:20, verbose=false);

1 thread:

21.393467 seconds (87.84 M allocations: 3.032 GiB, 3.15% gc time)
10.813218 seconds (52.30 M allocations: 1.065 GiB, 1.53% gc time)

4 threads

13.381170 seconds (88.03 M allocations: 3.043 GiB, 4.44% gc time)
3.201153 seconds (52.76 M allocations: 1.079 GiB, 6.05% gc time)

8 threads

12.632355 seconds (89.10 M allocations: 3.089 GiB, 4.85% gc time)
2.277163 seconds (52.36 M allocations: 1.069 GiB, 5.31% gc time)

[DanielEWeeks]

In your instructions, you state

Note: If one is running IHT on Matrix{Float64}, BLAS should NOT run with multiple threads (execute BLAS.set_num_threads(1) before running IHT). This prevents oversubscription.

but if I have only loaded MendelIHT, then the suggested command does not work:

julia> BLAS.set_num_threads(1)
ERROR: UndefVarError: BLAS not defined
Stacktrace:
 [1] top-level scope
   @ REPL[97]:1

I thought this meant that I needed to install the BLAS package, which surprised me as MendelIHT is running fine.

Looks like I just need to load LinearAlgebra before doing BLAS.set_num_threads(1), as this works fine:

julia> using LinearAlgebra
julia> BLAS.set_num_threads(1)

Consider clarifying this in your instructions re how to turn threading on.

Thanks,
Dan