The normals would be generated based on the gradient of the density function at the same time that you get the intersection points between the edges and the surface.  If it's something simple and closed-form like a sphere then you can calculate the normals analytically, but with noise you'll need to take samples.

You have the next steps in the wrong order.  First, you generate a vertex for each *cell* that exhibits a sign change.  The QEF you are minimizing is simply the total (squared) distance to each of the planes that are defined by the intersection point/normal pairs for that cell.  Then you walk through the edges that exhibit sign changes and create a quad using the four adjacent vertices (which are guaranteed to have been generated in the last step).

Now, my biggest hurdle in implementing this was solving the QEF.  I actually came up with a simple iterative solution that will run well on (for example) a GPU in parallel.  Basically, you start the vertex in the centre of the cell.  Then you average all the vectors taken from the vertex to each plane and move the vertex along that resultant, and repeat this step a fixed number of times.  I found moving it ~70% along the resultant would stabilize in the least amount of iterations.