How can I make regressions with "big data" faster in R? So I'm running a regression in R, with the following formula:
lm(y ~ x1 + x2 + factor(x3))
The issue is that x3 has 10000 levels, thus there are slightly over 10000 predictors. 
The regression runs fast when there are 30000 observations, but is excruciatingly slow when there are 300000 observations. Is there any way I can make it faster?
 A: Three ways, One is that the lm function is doing a LOT of stuff under the hood. It computes parameter values, stadnard deviations, p values, and many many more things. If you want to see this, call str(lm0) where lm0 if you linear model and you'll see a lot of things. So the first thing is to just write a much smaller function that just does what you want. Whether that be parameter values or values + pvalues. 
The second thing I would jump to is a different matrix factorization first. Linear regression in R can be sped up by doing a Gram-Schmidt facotrization of your matrix first before you regress. (https://cran.r-project.org/web/packages/matlib/vignettes/gramreg.html)
And finally, as for optimization algorithms, gradient descent can be faster. This is because R tries to find an analytic solution to your regression problem. But in larger data contexts, doing this is incredibly costly. So approximate solutions like gradient descent trade an exact solution that takes forever to compute with an approximate solution that is much quicker to solve. However you can do much better than gradient descent. Gradient Descent is a generic optimization algorithm that doesn't take advantage of the unique convex structure of the optimization problem. An algorithm that does is the Iteratively Reweighted Least Squares (IRLS) algorithm. This can be faster as long as you modify the glm.control parameters as the default parameters will be slower than the default lm function. And finally as someone mentioned in the comments above, the glmnet package also has highly efficient algorithms for fitting regression parameters.
Something that has also helped me quite a bit has been the Microsoft R (formerly  Revolution R) program. This implementation of R is said to be focused on bringing multithreded math libraries to R. And while knowledge of multithreading is outside of my expertise, I have gotten noticeable speed increases from this. The related Microsoft R Client also boasts better compatibility with spark and hadoop, which might be of help as well.  
A: A lot of fine comments and answers here, but is it all really necessary? This is a very simple model. Let's rewrite it as
$$
y = \beta_1(x_1 - \mu_1) + \beta_2(x_2 - \mu_2) + \phi_f + \epsilon
$$
where $\phi_f$ is the offset due to factor $f$. Clearly, $\phi_f$ is the average value of the response when the independent variables are at their mean values.
With the vast number of factor levels you have, it makes sense to treat $\phi_f$ as a random effect, so don't use glm or lm Using lmemight still cost a lot of time, but would you really lose that much by fitting the linear model and the factor model independently? If the sample is random and the number of items at each level is roughly the same, I'm not sure that you would lose much by estimating the random effects as the average over the $y$'s at each level. Then do the regression without the factors -- the factor effects will basically be rolled into the error $\epsilon$.
You won't get exactly the same answer as if you did it the right way, but realistically, is the extra computing time worth it, if that's causing you serious grief.
Another option, which I have never tried, but might help, would be to get the naive estimates I described above, and then do one or two iterations through the EM algorithm towards improved versions, but don't attempt to go the whole way to convergence with all that data and all those parameters. 
