Is there a closed-form solution for
$$\max_D \operatorname{Tr}(ADBD)$$
where $D$ is a $N\times N$ diagonal matrix with $m<N$ number of $1$'s and the rest are $0$'s, and $A$ and $B$ are real positive semi-definite matrices. $m$ is given. I am not interested in obtaining $\arg \max$, but only $\max$.
For any of your diagonal matrices $D$, let $J:=J_D$ be the set such that $D_{i,j}=1(i=j\in J)$ for all $i,j$ in the set $[N]:=\{1,\dots,N\}$, where, for any matrix $M$, its $(i,j)$-entry is denoted by $M_{i,j}$. Then $J$ is in $\binom{[N]}m$, the set of all subsets of $[N]$ of cardinality $m$.
So, for $J:=J_D$, $$\text{Tr}(ADBD) =\sum_{i,j,k,l\in[N]}A_{i,j}1(j=k\in J)B_{k,l}1(l=i\in J) \\ =\sum_{i,j\in J}A_{i,j}B_{j,i}=\text{Tr}(A_J B_J),$$ where $M_J:=[M_{i,j}\colon i,j\in J]$.
So, $\max_D\text{Tr}(ADBD)$ is the maximum of $\text{Tr}(A_J B_J)$ over all $J\in\binom{[N]}m$.
I can reduce this discrete optimization problem to a continuous optimization problem in many cases by using my own techniques, but the continuous optimization problem may return a higher value than the discrete optimization problem in some cases. Sometimes the continuous optimization problem finds a solution that does not exist in the discrete case.
The OP has observed that if R,S are matrices, then $\text{Tr}(RS)=\sum(R\odot S^T)$ and therefore $\text{Tr}(ADBD)=\sum(D(A\odot B^T)D)$.
The problem of maximizing $\sum(DCD)$ where $D$ is a 0-1 diagonal matrix with $\text{Tr}(D)=m$ is NP-complete since the clique problem in graph theory is a special case of this problem. But as with many NP-complete problems, one can develop heuristic algorithms for maximizing $\sum(DCD)$ such as hill climbing algorithms or evolutionary algorithms.
The expression $\max_D\sum(DCD)$ is a limit case of the expression $\max_D\rho(DCD)$ where $\rho$ denotes the spectral radius. If $N_n$ is the $n\times n$-matrix where each entry in $N_n$ is $1$, and $C$ is a real matrix, then $\sum(C)=\lim_{t\rightarrow 0}\frac{\rho(nt\cdot C-N_n)-n}{t}=\frac{\partial}{\partial t}\rho(nT\cdot C-N_n)|_{t=0}$ (this limit can be derived from the formula for the gradient of the eigenvalues of a matrix), and we have a similar expression for $\max_{D}\sum(DCD)$ in terms of the spectral radius.
Suppose that $K$ is either the field of real numbers or complex numbers (we will need to make modifications to generalize everything to quaternionic matrices) and $(A_1,\dots,A_r),(B_1,\dots,B_r)\in M_n(K)^r$. Then define the $L_2$-spectral radius of $(A_1,\dots,A_r)$ by $\rho_2(A_1,\dots,A_r)=\rho(A_1\oplus\overline{A_1}+\dots+A_r\oplus\overline{A_r})$. If $(A_1,\dots,A_r)\in M_m(K),(B_1,\dots,B_r)\in M_n(K)$, then define the $L_2$-spectral radius similarity between $(A_1,\dots,A_r)$ and $(B_1,\dots,B_r)$ is the quantity $$\|(A_1,\dots,A_r)\simeq(B_1,\dots,B_r)\|_2=\frac{\rho(A_1\oplus\overline{B_1}+\dots+A_r\oplus\overline{B_r})}{\rho_2(A_1,\dots,A_r)\rho_2(B_1,\dots,B_r)}.$$
Define the $L_{2,d}$-spectral radius of $A_1,\dots,A_r$ by setting $$\rho_{2,d}^K(A_1,\dots,A_r)=\max\frac{\rho(A_1\oplus\overline{X_1}+\dots+A_r\oplus \overline{X_r})}{\rho_2(X_1,\dots,X_r)}:X_1,\dots,X_r\in M_d(K)\}.$$
If $(X_1,\dots,X_r)\in M_d(K)^n$ locally maximizes the quantity $$\frac{\rho(A_1\oplus \overline{X_1}+\dots+A_r\oplus\overline{X_r})}{\rho_2(X_1,\dots,X_r)},$$ then we say that $(X_1,\dots,X_r)$ is an $L_2$-spectral radius dimensionality reduction (LSRDR) of $(A_1,\dots,A_r)$. We observe that $(X_1,\dots,X_r)$ is an LSRDR of $(A_1,\dots,A_r)$ precisely when $(X_1,\dots,X_r)$ maximizes the similarity $\|(A_1,\dots,A_r)\simeq(X_1,\dots,X_r)\|$.
Empirical observation: Suppose that $(X_1,\dots,X_r)\in M_d(K)^r$ is an LSRDR of $(A_1,\dots,A_r)$. Then we can typically find matrices $R,S$ and some constant $\lambda$ where $X_j=\lambda RA_jS$ for all $j$ and where $RS$ is the identity matrix, and $SR$ is a (typically non-orthogonal) projection. In this case, let $P=SR$. The projection matrix $P$ is often unique in the sense that if we compute the LSRDR again and compute another projection matrix $P_1$, then we would get $P=P_1.$ Observe that $\rho_2(X_1,\dots,X_r)=\rho_2(PA_1P,\dots,PA_rP)$ and $\rho(A_1\otimes\overline{X_1}+\dots+A_r\otimes\overline{X_r})=\rho(A_1\otimes\overline{PA_1P}+\dots+A_r\otimes\overline{PA_rP})$.
Given $A_1,\dots,A_r\in M_n(K),X_1,\dots,X_r\in M_d(K)$, define a superoperator $$\Gamma(A_1,\dots,A_r;X_1,\dots,X_r):M_{n,d}(K)\rightarrow M_{n,d}(K)$$ by letting $$\Gamma(A_1,\dots,A_r;X_1,\dots,X_r)(X)=\sum_{k=1}^rA_kXX_r^*,$$ and set $$\Phi(A_1,\dots,A_r)=\Gamma(A_1,\dots,A_r;A_1,\dots,A_r).$$
Suppose now that $C=(c_{i,j})_{i,j}$ is a matrix with non-negative real entries. Let $a=(a_{i,j})_{i,j}$ be a matrix with $|a_{i,j}|^2=c_{i,j}$ for all $i,j$. For each $i,j$, let $A_{i,j}$ be the matrix where each entry except for the $i,j$ entry is $0$, but where the $(i,j)$ entry of $A_{i,j}$ is $a_{i,j}$. Suppose now that $(X_{i,j})\in M_d(K)^{n\times n}$ is an LSRDR of $(A_{i,j})_{i,j}$. Then the projection matrix $P$ is often a diagonal projection matrix where the diagonal entries are either $0$ or $1$.
If $v$ is a vector, then let $\text{Diag}(v)$ denote the square matrix whose diagonal entries are $v$, and if $V$ is a matrix, then let $\text{Diag}(V)$ denote the vector whose entries are the diagonal entries in $V$.
$\Gamma((A_{i,j})_{i,j},(PA_{i,j}P)_{i,j})=\Phi((PA_{i,j}P)_{i,j})$, and $\Phi((PA_{i,j}P)_{i,j})(X)=\text{Diag}(C\cdot\text{Diag}(X))$ for all matrices $X$. This means that $$\rho(\Gamma((A_{i,j})_{i,j},(PA_{i,j}P)_{i,j}))=\rho(\Phi((PA_{i,j}P)_{i,j})=\rho(PCP),$$ so
$$\max\{\rho\big(\text{Diag}(D)\cdot C\cdot\text{Diag}(D)\big):D\in\{0,1\}^n,\sum D=d\}:=\rho_{2,d}((A_{i,j})_{i,j})^{1/2}$$ if everything works well.
I may add details to this answer about how we can write $$\max\{\rho\big(\text{Diag}(D)\cdot C\cdot\text{Diag}(D)\big):D\in\{0,1\}^n,\sum D=d\}$$ as a continuous optimization problem without needing to use the spectral radius.
You can find more information about LSRDRs at my cryptocurrency site circcashcore.org/blog (I use LSRDRs and similar machine learning algorithms to evaluate the cryptographic security of block ciphers).
Added 4/26/2024
We observe that for sufficiently large $\delta$, we have $$\max\{\sum\text{Diag}(d)(C+I\delta)\text{Diag}(e):d,e\in\{0,1\}^n,\sum d=\sum e=m\}$$ $$=\max\{\sum\text{Diag}(d)(C+I\delta)\text{Diag}(d):d\in\{0,1\}^n,\sum d=m\}$$ $$=m\delta+\max\{\sum\text{Diag}(d)C\text{Diag}(d):d\in\{0,1\}^n,\sum d=m\}.$$
Therefore, the problem of finding the maximum element in the set $$\{\sum\text{Diag}(d)C\text{Diag}(d):d\in\{0,1\}^n,\sum d=m\}$$ is a special case of the problem of finding the maximum element in the set $$\{\sum\text{Diag}(d)C\text{Diag}(e):d,e\in\{0,1\}^n,\sum d=\sum e=m\},$$ but I also can rewrite this maximization problem in terms of what I call CLSRDRs.
Suppose that $A=(a_{i,j})_{i,j}$ is an $m\times n$ matrix with non-negative real entries. Then for $1\leq i\leq m$ and $1\leq j\leq n$, let $A_i,B_{i,j},C_j$ be the $1\times m,m\times n,n\times 1$ matrices respectively where each of the matrices $A_i,B_{i,j},C_j$ has at most one non-zero entry. The $(1,i)$ entry of $A_i$ is $1$, the $(i,j)$ entry of $B_{i,j}$ is a value $c_{i,j}$ where $|c_{i,j}|^2=a_{i,j}$ and the $(j,1)$ entry of $C_j$ is $1.$ Suppose now that $1\leq p\leq m,1\leq q\leq n$.
Let $$V=\max\frac{\rho(\Gamma((A_i)_i;(X_i)_i)\Gamma((B_{i,j})_{i,j};(Y_{i,j})_{i,j})\Gamma((C_j)_j;(Z_j)_j))}{\rho(\Phi((X_i)_i)\Phi((Y_{i,j})_{i,j})\Phi((Z_i)_i))^{1/2}}$$ subject to the condition that $X_i\in M_{1,p}(K),Y_{i,j}\in M_{p,q}(K),Z_j\in M_{q,1}(K)$ where $K\in\{\mathbb{R},\mathbb{C},\mathbb{H}\}$. Then in general, we can find matrices $S_1,S_2,R_2,R_3$ where $S_1\in M_{m,p}(K),S_2\in M_{n,q}(K),R_2\in M_{p,m}(K),R_3\in M_{q,n}(K)$ and $R_2S_1,R_3S_2$ are the identity matrices, $S_1R_2,S_2R_3$ are diagonal projection matrices, and where we reach our maximum $V$ when $X_i=A_iS_1,Y_{i,j}=R_2B_{i,j}S_2,Z_j=R_3C_j$.
In this case, $V^{1/2}=\max\{\sum\text{Diag}(d)A\text{Diag}(e):d,e\in\{0,1\}^n,\sum d=p,\sum e=q\}$ in many but not all instances.
