ENH - slow ``gram_cd_solver`` when fitted on sparse dataset

[Badr-MOUFAD]

As investigated in #59 through timing and benchmarks, fitting a gram_cd_solver on a sparse dataset comes with a big overhead of computing the gram matrix, as opposed to the dense case.

Sparse case

Dense case

It seems that it has to do more with scipy sparse matrix multiplication. More details about that (explanation, code snippets).

It would be beneficial to have a more efficient way to compute the gram matrix in the sparse case as it would speed up drastically the solver.