The Gram-Schmidt process is locally smooth, i.e. there is a small neighborhood around every basis where you don't encounter the problem you mention and where the orthonromal basis from produced depends smoothly on the input data. The problem that you describe appears only when you consider it globally, as was noted in other answers and comments.

There are other algorithms for QR decomposition apart from GS based on Householder or Givens rotations and they have some advantages (numerical stability, paralellizability). The appropriate search term would be "thin QR decomposition" or "reduced QR decomposition". 

Please note that QR decompositions actually produce orthogonal basis with a rather strong property: first $l$-tuples of vectors $\{q_1, \ldots, q_l\}$ generate the same subspace as $\{b_1, \ldots, b_k\}$  for all $l\leq k$. 

The method proposed by Tom Goodwillie is sometimes called "symmetric orthogonalization" or "SVD based orthogonalization". If you consider SVD of $B = U\Sigma V^t$, then the orthogonal basis is given by rows of $UV^t$. This solution has a remarkable property that it is the [nearest "orthogonal" matrix][1] to $B$. You can find more under the name [polar decomposition][2].

Kind of midway between these two approaches is QR decomposition with [pivoting or rank revealing QR decomposition][3] which tackles the problem of $BB^t$ having small eigenvalues (which corresponds to the basis being close to a degenerate one).


  [1]: https://en.wikipedia.org/wiki/Singular_value_decomposition#Nearest_orthogonal_matrix
  [2]: https://en.wikipedia.org/wiki/Polar_decomposition#Numerical_determination_of_the_matrix_polar_decomposition
  [3]: https://en.wikipedia.org/wiki/QR_decomposition#Column_pivoting