Given matrix $[D]\in \Bbb R^{m \times n}$ of spectral data composed of $n$ mixtures of two pure component spectra $\vec p$ and $\vec q\in\Bbb R^{m}$, scaled by matrix $[C]\in\Bbb R^{2 \times n}$ consisting of the concentrations of the two components in each mixture. The $m$-dimension represents the spectrum x-axis index and the $n$-dimension represents the sample mixture index. We have: $$[D]_{m\times n} = \begin{bmatrix}\vec p \\\vec q \end{bmatrix}^{\top}_{m\times2}[C]_{2\times n}$$
Imagine a set of spectra $[D]$ that have been normalized by the concentration of the second component and that we know the first pure component spectrum $\vec p$. Thus we have: $$[D]_{m\times n} = \begin{bmatrix}\vec p \\\vec q \end{bmatrix}^{\top}_{m\times2}\begin{bmatrix}\,\vec c\\1\cdots 1\end{bmatrix}_{2\times n}$$How can we obtain the least-squares estimates of $\vec q$ and $\vec c$?
The notation $[D]: \mathbb R \to \mathbb R^{m \times n}$ confuses me, so I suppose you've meant simply $D \in \mathbb R^{n \times m}$. I also will treat $p, q \in \mathbb R^m$ and $c \in \mathbb R^n$ as column vectors.
Let's simplify the matrix notation: $$ D = \begin{bmatrix} p^\top\\ q^\top \end{bmatrix}^\top \begin{bmatrix} c^\top\\ e^\top \end{bmatrix} = \begin{bmatrix} p & q \end{bmatrix} \begin{bmatrix} c^\top\\ e^\top \end{bmatrix} = p c^\top + q e^\top. $$ Here $e \in \mathbb R^n$ denotes a vector of ones. It is now obvious that the right hand side is linear both in $c$ and $q$.
The least squares solution minimizes the residual's sum of squares that is $$ E = \sum_{i=1}^m \sum_{j=1}^n (d_{ij} - p_i c_j - q_i)^2 \to \min_{c_j, q_i}. $$ This is a quadratic problem and the solution can be obtained from the optimality conditions: $$ 0 = \frac{\partial E}{\partial c_j} = \sum_{i=1}^m 2 (d_{ij} - p_i c_j - q_i) (-p_i) = -\sum_{i=1}^m p_i (d_{ij} - p_i c_j - q_i) = -2p^\top (D - pc^\top - qe^\top)\\ 0 = \frac{\partial E}{\partial q_i} = -\sum_{j=1}^n 2 (d_{ij} - p_i c_j - q_i) = -2 (D - pc^\top - qe^\top) e = 0. $$ Transposing the first equation we have $$ (D^\top - cp^\top - eq^\top) p = 0. $$ Now separating unknowns from the knowns we obtain a system of linear equaions: $$ (cp^\top + e q^\top) p = D^\top p\\ (pc^\top + qe^\top) e = D e $$ Using property $a^\top b = b^\top a$ when the product is scalar and noting that $a^\top b$ as a scalar can be safely moved across the product we get $$ (p^\top p) c + (e p^\top) q = D^\top p\\ (pe^\top) c + (e^\top e) q = De $$ In matrix form this system can be written as $$ \begin{bmatrix} p^\top p I & ep^\top\\ pe^\top & n I \end{bmatrix} \begin{bmatrix} c\\q \end{bmatrix} = \begin{bmatrix} D^\top p\\De \end{bmatrix} $$ This system's matrix is singular, but consistent. It can be checked by multiplying with $[e^\top\; -p^\top]$ on the left. The solution is also not unique. If $c_0, q_0$ is a solution then $$ c = c_0 + \alpha e, \quad q = q_0 - \alpha p. $$ also would be. $c_0, q_0$ may be obtained using the pseudoinverse matrix.
Note that for arbitrary $\alpha$ the difference between $D$ and $pc^\top + qe^\top$ will be the same due to $$ pc^\top + qe^\top = pc_0^\top + \alpha pe^\top + q_0 e^\top - \alpha pe^\top = pc_0^\top + q_0 e^\top. $$