What is a good way to import interpolating functions of large size?

Question

I was generating a set of eigenfunctions from a differential equation using NDEigensystem. Because of the computer time cost is too high, I thought it'd be better if I exported the data (which means the interpolating functions) to a file so I can import it anytime later without extra calculation.

Here is my code;

V[r_] = -(1/r) - (1.0415223038416566` E^(-0.9990999998788636` r))/r; 
α = 2.;
Module[{shift = 10, d = 2000, n = 20, ev, evShifted}, 
  {ev, ef} = 
    NDEigensystem[
     {shift f[r] + V[r] f[r] - 1/α f''[r], 
      DirichletCondition[f[r] == 0, True]}, 
     f, {r, 0, d}, n, 
     Method -> 
       {"SpatialDiscretization" -> 
          {"FiniteElement", {"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}, 
        "Eigensystem" -> {"Arnoldi", MaxIterations -> Infinity}}];
  evShifted = ev - shift]
Export["C:\\Users\\ASUS\\Documents\\TestDATA_efHQ.dat", ef, "Table"];

The total size of my exported file "TestDATA_efHQ.dat" is up to 4GB. When I tried to import it back, it drained my 12GB memory instantly and cost me almost 3 hours before the kernel failed.

I did try ReadList once, but it generated a stream which I don't know how to convert to an InterpolatingFunction, so I cut it before it's finished.

If anyone knows how to speed up this process or has suggestions about exporting and importing InterpolatingFunction objects, I'd be very grateful to hear what they have to say. I haven't try to export my data into some different format, such as binary, so if that would make this quicker, please let me know.

BTW, while I was using this method to do some integrations and calculate some function values, I found that for eigenfunctions with higher resolution, the method has better convergence and also gives higher accuracy than other methods I have tried.