What are the software limitations in all possible subsets selection in regression?  If I have a dependent variable and $N$ predictor variables and wanted my stats software to examine all the possible models, there would be $2^N$ possible resulting equations. 
I am curious to find out what the limitations are with regard to $N$ for major/popular statistic software since as $N$ gets large there is a combinatorial explosion.
I've poked around the various web pages for packages but not been able to find this information. I would suspect a value of 10 - 20 for $N$?
If anyone knows (and has links) I would be grateful for this information.
Aside from R, Minitab, I can think of these packages SAS, SPPS, Stata, Matlab, Excel(?), any other packages I should consider?
 A: I was able to generate all possible subsets using 50 variables in SAS.  I do not believe there is any hard limitation other than memory and CPU speed.
Edit
I generated the 2 best models for N=1 to 50 variables for 5000 observations.
@levon9 - No, this ran in under 10 seconds. I generated 50 random variables from (0,1)
-Ralph Winters
A: As $N$ gets big, your ability to use maths becomes absolutely crucial.  "inefficient" mathematics will cost you at the PC.  The upper limit depends on what equation you are solving.  Avoiding matrix inverse or determinant calculations is a big advantage.
One way to help with increasing the limit is to use theorems for decomposing a large matrix inverse from into smaller matrix inverses.  This can often means the difference between feasible and not feasible.  But this involves some hard work, and often quite complicated mathematical manipulations!  But it is usually worth the time.  Do the maths or do the time!
Bayesian methods might be able to give an alternative way to get your result - might be quicker, which means your "upper limit" will increase (if only because it gives you two alternative ways of calculating the same answer - the smaller of two, will always be smaller than one of them!).
If you can calculate a regression coefficient without inverting a matrix, then you will probably save a lot of time.  This may be particularly useful in the Bayesian case, because "inside" a normal marginalisation integral, the $X^{T}X$ matrix does not need to be inverted, you just calculate a sum of squares.  Further, the determinant matrix will form part of the normalising constant.  So "in theory" you could use sampling techniques to numerically evaluate the integral (even though it has an analytic expression) which will be eons faster than trying to evaluate the "combinatorical explosion" of matrix inverses and determinants.  (it will still be a "combinatorical explosion" of numerical integrations, but this may be quicker).
This suggestion above is a bit of a "thought bubble" of mine.  I want to actually test it out, see if it's any good.  I think it would be (5,000 simulations + calculate exp(sum of squares) + calculate weighted average beta should be faster than matrix inversion for a big enough matrix.)
The cost is approximate rather than exact estimates.  There is nothing to stop you from using the same set of pseudo random numbers to numerically evaluate the integral, which will again, save you a great deal of time.
There is also nothing stopping you from using a combination of either technique.  Use exact when the matrices are small, use simulation when they are big.  This is because in this part of the analysis.  It is just different numerical techniques - just pick the technique which is quickest!
Of course this is all just a bit of "hand wavy" arguments, I don't exactly know the best software packages to use - and worse, trying to figure out which algorithms they actually use.
A: I suspect 30--60 is about the best you'll get. The standard approach is the leaps-and-bounds algorithm which doesn't require fitting every possible model. In $R$, the leaps package is one implementation.
The documentation for the regsubsets function in the leaps package states that it will handle up to 50 variables without complaining. It can be "forced" to do more than 50 by setting the appropriate boolean flag.
You might do a bit better with some parallelization technique, but the number of total models you can consider will (almost undoubtedly) only scale linearly with the number of CPU cores available to you. So, if 50 variables is the upper limit for a single core, and you have 1000 cores at your disposal, you could bump that to about 60 variables.
A: Just a caveat, but feature selection is a risky business, and the more features you have, the more degrees of freedom you have with which to optimise the feature selection criterion, and hence the greater the risk of over-fitting the feature selection criterion and in doing so obtain a model with poor generalisation ability.  It is possible that with an efficient algorithm and careful coding you can perform all subsets selection with a large number of features, that doesn't mean that it is a good idea to do it, especially if you have relatively few observations.  If you do use all subsets selection, it is vital to properly cross-validate the whole model fitting procedure (so that all-subset selection is performed independently in each fold of the cross-validation).  In practice, ridge regression with no feature selection often out-performs linear regression with feature selection (that advice is given in Millar's monograph on feature selection).
