Is there any closed solution or algorithm for solving the following optimization problem based on $Z$, $X$ as an input?
$D \leftarrow \operatorname{argmin}\|X − DZ\|^2_F$ s.t $\|d_k\|_2^2 = 1, \forall k \in \{1,...K\}.$
where $d_k$ is kth column of $D$ and $F$ stands for Frobenius norm.
$Z \in \mathbb{R}^{K\times n}$, $X \in \mathbb{R}^{m\times n}$, and $D \in \mathbb{R}^{m\times K}$.
$\def\p#1#2{\frac{\partial #1}{\partial #2}}\def\D{\operatorname{Diag}}\def\S{\operatorname{Sym}}$Here is an unconventional approach to your problem.
For algebraic convenience, use a colon to denote the trace, i.e. $$A:B = {\rm Tr}(A^TB)$$ This Frobenius product acquires several nice properties from the underlying trace function, e.g. $$\eqalign{ A:A &= \big\|A\big\|_F^2 \\ A:B &= B:A = B^T:A^T \\ CA:B &= A:C^TB = C:BA^T \\ }$$ The objective function is then $$\phi(D) = \big(DZ-X\big):\big(DZ-X\big)$$ Next, introduce an unconstrained variable $U$ and use it to construct the $D$ matrix such that the constraint is automatically satisfied
$$\eqalign{ W &=\D\big(U^TU\big)^{-1/2}\quad\implies\quad W^2=\D\big(U^TU\big)^{-1}\\ D &= UW \quad\implies\quad I=\D(D^TD)\iff\big\|d_k\big\|^2_F=1 \\ }$$ where the function $\D(X)$ retains the main diagonal of its argument but sets all other elements to zero.
It will be useful to have the differentials of these quantities for later use (note that diagonal matrices behave much like scalars under differentiation) $$\eqalign{ dW &= -\tfrac 12\,\D\big(U^TU\big)^{-3/2}\D(U^TdU+dU^TU) \\&= -W^3\D\big(\S(U^TdU)\big) \\ dD &= dU\,W + U\,dW \\&= dU\,W-UW^3\D\big(\S(U^TdU)\big) \\ }$$ where the function $\;\S(X)\doteq\tfrac 12(X^T+X)$
Now we can calculate the gradient $\left(\p{\phi}{U}\right)$, set it to zero, and solve for the optimal matrix. $$\eqalign{ d\phi &= (DZ-X):dD\,Z \\ &= (DZ-X)Z^T:dD \\ &= M:dD \\ &= M:dU\,W - M:UW^3\D\big(\S(U^TdU)\big) \\ &= MW:dU - \D(W^3U^TM):U^TdU \\ &= \Big(M - UW^2\D(U^TM)\Big)W:dU \\ \p{\phi}{U} &= \Big(M - UW^2\D(U^TM)\Big)W \;\doteq\; 0 \\ M &= UW^2\D(U^TM) \\ U &= M\Big(W^2\D\big(U^TM\big)\Big)^{-1} \\ }$$ $\big[\,$The fact that $W$ is diagonal (and symmetric) was used to simplify several of the steps.$\,\big]$
While this isn't a closed form solution, it does provide the basis for an iterative solution.
After initializing $\,k=0,\;U_0=$ (a random matrix), perform the following sequence of calculations until convergence is achieved $$\eqalign{ W_{k+1} &= \D(U^T_kU_k)^{-1/2} \\ M_{k+1} &= U_kW_{k+1}ZZ^T-XZ^T \\ U_{k+1} &= M_{k+1}\Big(W_{k+1}^2\D\big(U^T_kM_{k+1}\big)\Big)^{-1} \\ k &= k+1 \\ }$$ The computation is very inexpensive since it's easy to invert (and extract roots of) a diagonal matrix. The desired solution is then $$D=\arg\min(\phi) = \lim_{k\to\infty}U_k W_k$$