In a linear regression problem with sparsity constraint, $P = (P_1, \cdots, P_N)^{T}$ is the column vector of the outputs, and $D = (d_{j, k})$ is the $(N \times M)$- dimensional matrix of inputs. The objective function is
$$\text{argmin}_{c \in \Bbb R^{M}} (\Vert P - Dc \Vert_2^2 + \lambda \Vert c \Vert_0)$$
in which $\Vert c \Vert_0 = \# \{j: c_j \neq 0\}$
I learnt that this problem is NP-hard, but I don't understand why.
What I think:
There are in total $N$ cases that we need to consider
$$N = C(M, 1) + C(M, 2) + \cdots + C(M, M)$$ in which $M$ is the dimension of the coefficient vector $c$. $C(M, n) = \begin{pmatrix} M\\ n \end{pmatrix} = \frac{M!}{(M-n)!n!} $.
For each tuple of selected features, we perform OLS (without taking into account the regularizer) and record the loss: $$L = RMSE + \lambda \Vert c \Vert_0$$
After doing $N$ such calculations, we can choose the tuple of features that yield the smallest $L$.
However, I don't know whether the result that obtained really is the solution to the original objective function since we don't take into account the effect of regularizer in the first place.
Or this algorithm does not make sense, instead of comparing $L = RMSE + \lambda \Vert c \Vert_0$ between tuples of different sizes, we can only compare the $RMSE$ of the set of feature tuples of the same size.