Let C={1,…,k} and N={k+1,…,m} be respectively the set of causal and non-causal SNPs. Consider the quadratic form Q=zˆT(R+RΣγR)−1zˆ=zˆT(Im+ΣγR)−1a inside the exponential function in N(zˆ|0,R+RΣγR), where a=R−1zˆ can be precomputed. We solve the linear system (Im+ΣγR)b=a for b by observing that the m – k elements in b corresponding to non-causal SNPs (γℓ=0) are bℓ=aℓ and the remaining elements result from solving a system of k equations (Ik+nsλ2RCC)bC=aC−nsλ2RCNaN, where RCC is the k × k correlation matrix of the causal SNPs and RCN the k×(m−k) submatrix of R corresponding to the correlations between the causal and non-causal SNPs. In addition, we observe that det(Im+ΣγR) is simply det(Ik+nsλ2RCC) after expanding with respect to the rows corresponding to non-causal SNPs. Computationally, these computations require one Cholesky decomposition with complexity O(k3) and provide thus a considerable saving compared to the naive way of decomposing the whole m × m matrix with complexity O(m3).
