ultimately the matrix you create in step 1 is where you're going to construct the global stiffness matrix, which is a direct summation of all the element stiffness matrices.
This global matrix is going to be used to solve a linear system $A x = b$, where $A$ is the global stiffness matrix you just constructed by "adding" all the element stiffness matrices in step 2 together.
You could initialize a dense matrix (i.e. have memory allocated for every possible $i,j$ combination of global indices) and use an implicit node numbering mapping function to satisfy the requirements of step 1, however many FEM problems are quite sparse (i.e. the number of entries in the matrix which are non-zero is small compared to the total number of entries), thus a sparse matrix which tries to allocate memory only for entries which are non-zero (along with some additional information which can be used to determine which entry has a given value).
To efficiently construct sparse matrices you typically want to pre-allocate all the memory it will need before-hand (similar to how you would pre-allocate a dense matrix and then just set the values to what you want). Doing so requires you to have some estimate of which entries in your matrix could ever possibly be non-zero, since it's better to slightly over-allocate than to under-allocate.
You can determine this by figuring out which $i,j$ global index combinations are a part of any local element stiffness matrix, and simply mark those entries as potentially non-zero, and everything else is guaranteed to be zero.
Then you can construct a sparse matrix using this information.
You don't have to take this exact approach necessarily, but this is one common way which finite element methods are implemented in practice because it's relatively easy to understand and efficient, and there are of course additional extensions/modifications you can make to this method.