I am wondering if it would be viable to transfer the SH coefficients calculation to a compute shader instead of doing it on CPU which for our engine requires a readback of the cube map texture. I am not entirely sure how to go about this since it will be hard to parallelize as each thread will be writing to all the coefficients. A lame implementation would be to have one thread running the entire shader but I think that's going into TDR territory.

Currently, I am generating an irradiance map but I am planning on switching to storing it inside a spherical harmonics because of the smaller footprint.

Does anyone have any ideas on how we can move this to the GPU or its just not a viable option?

Yes, you can absolutely do this. Projection onto SH is an integral of your signal multiplied with the SH basis functions, which means you can formulate it as a parallel reduction on the GPU. I actually implemented this many years ago for a lightmap baker, where the final lightmap coefficients were generated by rendering a hemicube per-texel and then projecting that onto SH. This is pretty much exactly what you're doing, except that you have a full cubemap and not just the upper half of one. That was one of the first compute shaders that I ever wrote, so I'm sure it's nowhere near optimal. It probably would have been more cache efficient to do the reduction in 2D tiles instead along rows, and I'm sure the shared memory access pattern was full of bank conflicts. But If you want to have a look it should help get you started on your own implementation.