Bayesian formulation of best subset regression We know Ridge is equivalent to using a Gaussian prior and Lasso is equivalent to using a double exponential prior. 
What is the Bayesian interpretation (implied prior) for the best subset regression? Or is it the case that there is no Bayesian formulation for this?
 A: This hinges on how you measure the "best" in the best-subset method (i.e., what metric you are using to compare the different models).  Most of the best-subset methods involve minimising some error metric composed of some negative multiple of the log-likelihood and a "penalty" term that may depend on the number of observations and the number of parameters in the model under consideration.  Bear in mind that the best-subset method chooses a model based on the (penalised) maximum likelihood estimator under the model, so it will give you a "best model" along with the corresponding MLE of the parameters of that model.

General form of the best-subset method: For example, suppose we are considering a model $\mathscr{M}$ with log-likelihood $\ell$ that depends on a parameter vector $\boldsymbol{\theta} \in \boldsymbol{\Theta}_\mathscr{M}$ with length $k$.  Suppose we observe the vector $\mathbf{x}$ composed of $n$ observations.  The error metric used in best-subset method will usually be of the form:
$$\text{Error}(\mathscr{M}) 
= \lambda (n,k) - \eta \max_\boldsymbol{\theta \in \boldsymbol{\Theta}_\mathscr{M}} \ell_\mathbf{x}(\boldsymbol{\theta})
= \min_\boldsymbol{\theta \in \boldsymbol{\Theta}_\mathscr{M}} [\lambda (n,k) - \eta \ell_\mathbf{x}(\boldsymbol{\theta})],$$
where $\lambda$ is a positive penalty function and $\eta>0$ is a positive multiplier of the maximised log-likelihood under the model.  The best-subset method chooses the model $\mathscr{M}^*$ with MLE $\boldsymbol{\theta}^*$ that minimises this error metric.  Thus, if we have some class $\mathscr{G}$ containing models, then we choose the model that satisfies:
$$\text{Error}(\mathscr{M}^*)
= \min_{\mathscr{M} \in \mathscr{G}} \text{Error} (\mathscr{M})
= \min_{\mathscr{M} \in \mathscr{G}} \min_\boldsymbol{\theta \in \boldsymbol{\Theta}_\mathscr{M}} [\lambda (n,k) - \eta \ell_\mathbf{x}(\boldsymbol{\theta})].$$

Bayesian posterior equivalence: The above method is an estimation method based on minimising an objective function, so it is worth investigating whether we can replicate it using maximisation of a posterior density under Bayesian analysis.  To do this, we find the appropriate form of the prior that gives us the equivalent minimisation, and then we check that this prior is a valid density.  If we have some prior $\pi(\mathscr{M},\boldsymbol{\theta})$ on the model and parameter then this leads to the corresponding posterior:
$$\pi(\mathscr{M},\boldsymbol{\theta}|\mathbf{x}) = L_\mathbf{x}(\boldsymbol{\theta}) \pi(\mathscr{M}, \boldsymbol{\theta}).$$
Now, if we set $\pi(\mathscr{M}, \boldsymbol{\theta}) \equiv \exp( - \lambda (n,k)/\eta)$ then the posterior maximum is:
$$\begin{aligned}
\max_\mathscr{M, \boldsymbol{\theta}} \pi(\mathscr{M},\boldsymbol{\theta}|\mathbf{x}) 
&= \max_\mathscr{M, \boldsymbol{\theta}} L_\mathbf{x}(\boldsymbol{\theta}) \pi(\mathscr{M}, \boldsymbol{\theta}) \\[6pt]
&= \max_\mathscr{M \in \mathscr{G}} \max_\boldsymbol{\theta \in \boldsymbol{\Theta}_\mathscr{M}} L_\mathbf{x}(\boldsymbol{\theta}) \pi(\mathscr{M}, \boldsymbol{\theta}) \\[6pt]
&= \max_\mathscr{M \in \mathscr{G}} \max_\boldsymbol{\theta \in \boldsymbol{\Theta}_\mathscr{M}} [\log \pi(\mathscr{M}, \boldsymbol{\theta}) + 
\eta \ell_\mathbf{x}(\boldsymbol{\theta})] \\[6pt]
&= \min_\mathscr{M \in \mathscr{G}} \min_\boldsymbol{\theta \in \boldsymbol{\Theta}_\mathscr{M}} [- \eta \log \pi(\mathscr{M}, \boldsymbol{\theta}) - \eta \ell_\mathbf{x}(\boldsymbol{\theta}) ] \\[6pt]
&= \min_\mathscr{M \in \mathscr{G}} \min_\boldsymbol{\theta \in \boldsymbol{\Theta}_\mathscr{M}} [\lambda (n,k) - \eta \ell_\mathbf{x}(\boldsymbol{\theta}) ] \\[6pt]
&= \text{Error}(\mathscr{M}^*). \\[6pt]
\end{aligned}$$
Thus, we can see that the best-subset method is equivalent to the maximum a posteriori (MAP) estimator using the prior:
$$\pi(\mathscr{M}, \boldsymbol{\theta}) \equiv \exp \bigg( - \frac{\lambda (n,k)}{\eta} \bigg).$$
Now, obviously this equivalence is only going to be valid if this function is indeed a valid probability density function over the class of models and parameters (i.e., it must sum to one and it should not depend on $n$).  This imposes some strict requirements on the penalty function $\lambda$, which in general can depend on $n$ and $k$.  Since $n$ depends on the observed data, if the function depends on this value then we have information from the data in the prior and so this is not a strict Bayesian analysis.  Moreover, if this prior does not sum to one then it is not a valid density and so the equivalence does not hold.  In this case, the only way we can obtain a Bayesian equivalent is to move some of the "prior" weight into the likelihood function, and this means that the equivalent Bayesian model uses a different likelihood function to the best-subset method.
In some cases, such as when using the best-subset method using AIC, the above "prior form" does not depend on $n$ but it also doesn't generally sum to one (i.e., it is not a valid density).  In this case it is possible to alter the Bayesian analysis by taking a scaling constant that depends on $k$ out of the prior (to make it sum to one) and putting it into the likelihood function.  Since $k$ depends on the parameter vector, this alters the likelihood function, and so it no longer corresponds to the likelihood under the best-subset method.  Nevertheless, you obtain an "equivalence" of sorts, using likelihood functions that differ by a scaling value that depends on the length of the parameter vector.
