Will this work ?

1. Store the adjacency list of each node in a set
2. Start from rhe node with the least degree to the most degree (denote this node as $$$u$$$) and do the following :
3. Iterate $$$v$$$ from $$$1$$$ to $$$N$$$. While edge $$$(u,v)$$$ is not added, add it to the graph (to check it, the complexity is $$$O(\text{log }N)$$$

Agree this approach will produce a less random. But how would you define "choice of edges should be provably random"?

You can generate $$$1$$$ edge at a time. To do that generate a random number between $$$1$$$ and $$$\frac{n \cdot (n - 1)}{2} - m$$$, get an edge by index with a binary search and then update everything with some data structures. One operation can be maintained in $$$O(\log n)$$$, however it has a bad constant factor.

Edit: actually, I can only see how to do in $$$O(\log ^2 n)$$$ because we need to get $$$k$$$-th element that is absent in some adjacency list, and I know how to dynamically do that only with ordered sets. There is probably a much better way.

You can achieve $$$O(\log n)$$$ complexity with that method if you use dynamic (implicit) segment trees instead of ordered sets, however I have a feeling that it would work even slower, and it uses $$$O(m \log n)$$$ memory.

Compiling what others have said, this is I think the best you can expect:

1. Generate a random sorted subset of size $$$K$$$ of numbers from $$$1$$$ to $$$N^2-M$$$ in $$$O(K \log K)$$$ with low constant.
2. Do two pointers to skip the $$$M$$$ edges while adjusting in linear time.

I've also included some code, as I think in this case it's more suggestive. Please note that I haven't tested the code whatsoever and there might be bugs, but it highlights the general idea.

vector<pair<int, int>> Generate(
    int n, int k, 
    vector<pair<int, int>>& edges, // expect all edges (a, b) to have b < a 
    mt19937_64& rng
) {
  // Generate k random indices from 1 to n^2 - m
  // We will afterwards sort them and make sure that they're distinct after
  // sorting by adding the number of smaller elements.
  ll sz = 1LL * n * (n - 1) / 2 - edges.size();
  vector<ll> indices(k);
  for (int i = 0; i < k; ++i)
    indices[i] = rng() % (sz--);
  sort(indices.begin(), indices.end());
  sort(edges.begin(), edges.end());
  
  // We skip existing edges via two pointers technique.
  // This part might be tricky to understand; it's best to think about:
  // 1. We leverage the bijection (a, b) <-> idx by the relationship
  //    idx = a * (a - 1) / 2 + b
  // 2. The `edges` array comes in (a, b) form, and the `indices` 
  //    array comes in idx form
  // Once you understand that, the "weird code" is just converting between
  // the two representations on the fly.
  vector<pair<int, int>> result;
  int ptr = 0, bias = 0, ea = 0;
  for (auto idx : indices) {
    idx += (bias++);
    while (ptr < (int)edges.size()) {
      // Check if (a, b) <= idx; if so, increase bias, otherwise break.
      auto [a, b] = edges[ptr];
      assert(b < a);
      ll jdx = 1LL * a * (a - 1) / 2 + b;
      if (jdx <= idx) ++ptr, ++bias, ++idx;
      else break;
    }
    // (ea, eb) is the edge to add now.
    while (1LL * (ea + 1) * ea / 2 <= idx) ++ea;
    int eb = idx - ea;
    assert(0 <= eb && eb < ea);
    result.emplace_back(ea, eb);
  }

  return result;
}

https://gist.github.com/theabbie/0be86edab481af581c450a4cfeef6d69

Implemented this random graph generator based on ideas suggested.