# Least Squares derivation for Multivariate Multiple Regression

I have a multivariate multiple regression problem of the form:

$$Y = X \beta$$

Where $$Y$$ is an $$n \times m$$ matrix, $$X$$ is $$n \times p$$, and $$\beta$$ is $$p \times n$$.

I understand the derivation of the least squares solution to the standard regression problem, where $$Y$$ and $$\beta$$ are vectors rather than matrices, but am struggling to derive the solution for this multiple case. Specifically, I'm not sure how to translate the loss function to this higher-dimensional case.

All the sources I've read so far only provide derivations for the standard multivariate regression problem. This makes me wonder if the solution to the multivariate multiple case is to split the problem into one response at a time, then aggregate the solutions. Is this the actual approach used by computational solutions (such as lsfit in R)? And if not, what is the derivation of the least squares solution?