This is the 0-1 quadratic knapsack problem, which is NP-hard. The binary decision variables $x_i$ indicate whether $i\in S$, and the knapsack capacity is $k$.
You can solve it via integer linear programming as follows. For $i<j$, let binary decision variable $y_{i,j}$ represent $x_i x_j$. The problem is to maximize $\sum_i M_{i,i} x_i + \sum_{i<j} (M_{i,j}+M_{j,i}) y_{i,j}$ subject to \begin{align} y_{i,j} &\le x_i &\text{for $i<j$} \tag1 \\ y_{i,j} &\le x_j &\text{for $i<j$} \tag2 \\ y_{i,j} &\ge x_i + x_j - 1 &\text{for $i<j$} \tag3 \\ \sum_i x_i &\le k \tag4 \end{align} Constraint $(1)$ enforces $y_{i,j} \implies x_i$. Constraint $(2)$ enforces $y_{i,j} \implies x_j$. Constraint $(3)$ enforces $(x_i \land x_j) \implies y_{i,j}$. Constraint $(4)$ enforces $|S| \le k$.
If $M_{i,j} \ge 0$, you can omit $(3)$, which will naturally be satisfied because of the objective.