In principle, if the best subset can be found (i.e. for large enough datasets or if the solution is extremely sparse or if there are nonnegativity constraints which can help to identify the best subset), it is indeed better than the LASSO, in terms of (1) selecting the variables that actually contribute to the fit, (2) not selecting the variables that do not contribute to the fit, (3) prediction accuracy and (4) producing essentially unbiased estimates for the selected variables. One recent paper that argued for the superior quality of best subset over LASSO is that by Bertsimas et al (2016) "Best subset selection via a modern optimization lens". Another older one giving a concrete example (on the deconvolution of spike trains) where best subset was better than LASSO or ridge is that by de Rooi & Eilers (2011). (Well that one actually approximates best subset using an iterative adaptive ridge regression approach as in Frommlet & Nuel (2016) and as used in my L0glm package that is in development)
The reason that the LASSO is still preferred in practice is mostly due to it being computationally much easier to calculate and the coefficient paths being quite stable. Best subset selection, i.e. using an $L_0$ pseudonorm penalty, is often seen as a combinatorial problem, and is then NP hard, whereas the LASSO solution is convex and easy to calculate over a regularization path using pathwise coordinate descent. In fact, the LASSO ($L_1$ norm penalized regression) is the tightest convex relaxation of $L_0$ pseudonorm penalized regression / best subset selection and so is uniquely defined (bridge regression, i.e. $L_q$ norm penalized regression with q close to 0 would in principle be closer to best subset selection than LASSO, but this is no longer a convex optimization problem, and so is quite tricky to fit).
To reduce the bias of the LASSO one can use derived multistep approaches, such as the adaptive LASSO (where coefficients are differentially penalized based on a prior estimate from a least squares or ridge regression fit, e.g. just using the LASSO selected variables and then using penalty weights = 1/abs(LASSO coefficients of 1st LASSO fit)) or relaxed LASSO (a simple solution being to do a least squares fit of the variables selected by the LASSO). In comparison to best subset, LASSO tends to select slightly too many variables though. Best subset selection is better, but harder to fit.
That being said, there are also efficient computational methods now to approximate best subset selection / $L_0$ penalized regression in a very stable way. One is using the adaptive ridge approach described in the paper
"An Adaptive Ridge Procedure for L0 Regularization" by Frommlet & Nuel (2016) (this approach has a nice empirical Bayesian interpretation), as implemented in the l0ara package (there referred to as broken adaptive ridge) and my own L0glm package, still in development. Here one is approximating the non-convex best subset objective by a series of convex optimization problems, which then still has a unique solution at convergence. Another one is the approach for fast best subset selection as described in the abess package, which executes in polynomial time. And there is also the L0Learn package for sparse learning & approximating the L0 penalty, also potentially in combination with L2 penalty regularisation to deal better with collinearity. All are well worth checking out. Traditional combinatorial approaches to approximate best subset, such as used in the glmulti, leaps or bestglm packages I would not recommend - they do not scale to large problem sizes & tend to be unstable.
Note that also for best subset selection you'll still have to use either cross validation or some information criterion to determine what number of predictors gives you either the best prediction performance / explanatory power for the number of variables in your model (e.g. based on the AIC or for high dimensional settings the GIC criterion) or approximates best the true support (e.g. based on BIC or for high dimensional settings eBIC or MBIC), which is essential to avoid overfitting.
The paper "Extended Comparisons of Best Subset Selection, Forward Stepwise Selection, and the Lasso" by Hastie et al (2017) provides an extensive comparison of best subset, LASSO and some LASSO variants like the relaxed LASSO, and they claim that the relaxed LASSO was the one that produced the highest model prediction accuracy under the widest range of circumstances, i.e. they came to a different conclusion than Bertsimas, but that was also because they mainly looked at predictive performance and not at consistency of variable selection. So the conclusion about which is best depends a lot on what your aim is, e.g. achieving highest prediction accuracy, or achieving variable selection consistency and being best at picking out relevant variables and not including irrelevant ones or recovering the true support size most accurately. Ridge regression e.g. typically selects way too many variables but the prediction accuracy for cases with highly collinear variables can nevertheless be really good. If the true solution is really sparse there are cases where the optimal level of regularisation that results in optimal predictive performance and optimal support recovery coincide though - that is of course ideal. But of course it also depends how you would like to choose the tradeoff between Type I & Type II errors. Ideally you would like to have both a low number of false positives and a low number of false negatives, but that may only be possible with relatively large datasets & when the true solution is very sparse.
For a problem with many more cases than variables like you describe it is plain clear best subset selection is the preferred option though.
In my tests I found that an iterative adaptive ridge approach / broken adaptive ridge / L0glm approach has a series of advantages: (1) it nicely approximates best subset, (2) does so in a stable way (i.e. it has stable coefficient paths over the level of regularisation), (3) has good variable selection consistency (e.g. the union of selected variables over bootstrapped datasets typically converges quickly to a stable set, whereas LASSO would keep on adding spurious variables, which could only be mitigated by selecting the variables occurring in more than 80% of the bootstrap replicates, as in the bolasso method), (4) deals well with collinearity and high dimensional problems (e.g. with several perfectly collinear variables in the model their effect would just be spread equally - it does grouped variable selection like elastic net), (5) resolves common problems one sometimes runs into with unregularised approaches (e.g. complete separation in logistic regression), (6) has the advantage that only variables with small or zero effect are removed and that variables with large effect will stay in the model and will be unbiased, implying that omitted variable bias will be small, (7) is fast and scales well to large problems with millions of variables (computational speed is similar to glmnet), (8) allows one to choose the optimal level of regularisation a priori based on what information criterion one would like to optimise, (9) allows one to do inference via bootstrapping for both low & high dimensional problems and (10) allows one to specify prior beliefs based on the starting values given (which correspond to the square root of the SD of the Gaussian prior implied in the first iteration of the algorithm) and the prior mean (normally zero, but can also be specified to be nonzero for particular parameters if nonzero centered ridge regression updates are used). Only in high dimensional settings with an extremely sparse true solution I found MCP & SCAD penalized regression to perform better (these can be fit using ncvreg or ordinis), but even that could be covered by using MCP or SCAD penalized regression estimates as warm start starting values in an iterative adaptive ridge L0glm approach, in which case L0glm can find an even better solution still & remove some false positives.
But all this is talking about particular approximations of best subset that are uniquely determined, and are calculated either using a nonconvex penalty (in MCP and SCAD) or as a sequence of convex optimization problems (in iterative / broken adaptive ridge / L0glm). Many best subset algorithms approach the problem as a combinatorial optimization problem, and they suffer from poor stability. That, combined with the high computational cost, is probably the reason many stay away from them. But given that in my tests iterative / broken adaptive ridge / L0glm models were superior to LASSO in terms of predictive performance, biasedness of the coefficients and support recovery, stable best subset approximations should imho receive a lot more attention. See e.g. here for some examples and benchmarks that I made relating to the single or multichannel deconvolution of a blurred sparse spike train, where I compared some different approaches (L0glm, abess, L0Learn, glmnet & ordinis) (this is a problem where for single channel n=p, so in that sense a relatively small dataset, though with nonnegativity constraints on the coefficients), or a high dimensional example (with n=500 and p=1000000 variables and 50 contributing variables, so with a very sparse solution and also with nonnegativity constraints) where nonnegative MCP penalized regression fit using ordinis or L0glm with MCP as initialisation performed very well & where in both cases these approximations of best subset performed better than LASSO, both in terms of predictive performance and in terms of support recovery / variable selection consistency (ending up with just a couple of false negatives & 1 false positive), whilst also scaling well to large problems. Both of those examples are ones with nonnegativity constraints on the coefficients, which is a situation that makes recovery of the best subset a lot easier, as it acts as an effective regularisation method in itself - in fact, under that circumstance a hard thresholded version of the nonnegative LASSO outperforms plain nonnegative LASSO in support recovery and unconstrained LASSO in predictive performance, and even a hard thresholded version of nonnegative least squares would then outperform the plain LASSO.
