Suppose we have some $n \times m$ matrix $R$ and we want to find non-negative decomposition on matrices $P$ of dimension $n \times d$ and $d \times m$-matrix $Q$. But since exact decomposition usually does not exist we have to minimize functional $L$. Suppose we fixed non-negative matrix $Q$.
$$ L(R, PQ) \rightarrow \min_{P \geq 0} $$
and
$$ L(R, PQ) = \|R - PQ\|^2 = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \left( r_{ij} - \sum_{k = 1}^{d} p_{ik} q_{kj} \right)^2. $$
Find explicit solution of this problem.
I think to use Kuhn–Tucker theorem but we studied it only in continuous form (for integrals), so I little bit confused.
I am reasonably sure there is no closed form solution to this problem. But this is, in fact, just a nonnegative least squares problem (NNLS). It may not immediately look like it, but thanks to a little Kronecker product magic we can rewrite the objective function in this fashion: $$\|R-PQ\|_F^2 = \|\mathop{\textrm{vec}}(R)-(Q^T\otimes I)\mathop{\textrm{vec}}(P)\|_2^2$$ where $\text{vec}$ is the matrix-to-vector conversion operation (e.g., the columns stacked one upon another), and $\otimes$ is the Kronecker product. So any software that can solve NNLS problems can handle this, now that it is converted to vector form. If $P$ were not constrained negatively, then $P=RQ^\dagger$, where $Q^\dagger$ is the Moore-Penrose pseudoinverse of $Q$. If $QQ^T$ is full rank, this reduces to $P=RQ^T(QQ^T)^{-1}$.
Once you see it as just nonnegative least squares then you have a wealth of information available to you. It can be solved using any quadratic programming engine, for instance, or with specialized methods for NNLS like the active set method offered in the Wikipedia article. The Kronecker structure of the coefficient matrix $Q^T\otimes I$ can be used to accelerate many of the methods you might choose.
In fact, if you look at that NNLS Wikipedia page, you will see that non-negative matrix factorization is one of the cited applications.