You are artificially sparsifying the rhs. Your function grad_sparse already returns the rhs as is. You could then even save the additional memory allocation due to gradinds, if you left the rhs dense. Also, you would save the allocation of the index vector for the result, and, if you don’t need grad_vec, you could even solve in-place!
One challenge with the backslash is that it first determines an appropriate factorization for hess_mat. Then, depending on whether it is triangular or even diagonal or neither, it uses different “factorizations”. Unfortunately, this may yield either dense or sparse return values, as you noticed in the github issue.