this seems some unnecessary operation to me. What about (python)
        # gp book 3.27, matrix inversion lemma on
        # (K^-1 +W)^-1 = K -KW^0.5 B^-1 W^0.5 K
        C = (K.T * w_sqrt).T
        C = solve_triangular(L, C, lower=True)
        C = solve_triangular(L.T, C, lower=False)
        C = (C.T * w_sqrt).T
        C = K.dot(C)
        C = K - C