Minimizing the Frobenius Norm I would like to minimize the following expression with respect to matrix X:
$$\left \| A-BX \right \|_{F}$$
where A and B matrices are given and all the matrices have positive integer elements.
Any hints?
 A: Ah, I didn't realize that $X$ was required to be positive as well.
Are you sure a clean solution exists?
Since $BX$ is monotonic in all entries of $X$, dynamic programming can solve the problem in time pseudopolynomial in the magnitudes of the entries of $A$; this approach will only be practical for small $A$ however. 
Branch and bound codes like MIQPBB will work, and would be the first thing I'd try if I had to solve the problem in practice. See also this paper for an algorithm that alleges to take advantage of the convexity of the problem.
Does the nonnegativity lead to even better specialized algorithms? It's possible, though I don't see it at the moment. Your problem is closely related to other notable optimization problems, such as a modified closest vector problem: given a point in the first orthant and $n$ vectors in that orthant, what positive linear combination of those vectors brings you closest to the point? I did a quick literature search and didn't turn up much of relevance, but it could be out there...
A: The specific problem you have appears to be trivial for state-of-the-art MIQP solvers (i.e., branch and bound/cut). It is solved to 5-6 digits of accuracy in 0s. There might be smarter ways to solve the problem, but I reckon not much smarter, as it is a combinatorial problem by nature. Heuristics might work very well though, in case you don't care about optimality certificates.
Your data distribution was not given, but based on what you have told us, I've implemented a test with the MATLAB Toolbox YALMIP below (disclaimer, developed by me). The MIQP solvers I tried were Gurobi and Cplex, and the built-in solver bnb (which is absolutely terrible, but works ok here)
T = 10;
A = ceil(rand(8,3)*1e4);
B = ceil(randn(8,T)*100);
X = intvar(T,3,'full');
Z = A-B*X;
Constraints = [X>=0];
Objective = Z(:)'*Z(:);
optimize(Constraints,Objective)

BTW, if you have large numbers in both $A$ and $B$, I would work with a model which avoids squaring and multiplying that data
Z = intvar(8,3,'full');
Constraints = [Z == A - B*X, X>=0];
Objective = Z(:)'*Z(:);
optimize(Constraints,Objective)

Of course, one could use the original form
Constraints = [X>=0];    
Objective = norm(A-B*X,'fro')
optimize(Constraints,Objective)

but that leads to a MISOCP, which is solved much slower by the solvers above.
